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Abstract.  Approximate  separable  representations  of  Green’s  functions  for  differential 
operators  is  a  basic  and  an  important  aspect  in  the  analysis  of  differential  equations  and  in 
the  development  of  efficient  numerical  algorithms  for  solving  them.  Being  able  to  approx¬ 
imate  a  Green’s  function  as  a  sum  with  few  separable  terms  is  equivalent  to  the  existence 
of  low  rank  approximation  of  corresponding  discretized  system.  This  property  can  be 
explored  for  matrix  compression  and  efficient  numerical  algorithms.  Green’s  functions  for 
coercive  elliptic  differential  operators  in  divergence  form  have  been  shown  to  be  highly 
separable  and  low  rank  approximation  for  their  discretized  systems  has  been  utilized  to 
develop  efficient  numerical  algorithms.  The  case  of  Helmholtz  equation  in  the  high  fre¬ 
quency  limit  is  more  challenging  both  mathematically  and  numerically.  In  this  work,  we 
develop  a  new  approach  to  study  approximate  separability  for  the  Green’s  function  of 
Helmholtz  equation  in  the  high  frequency  limit  based  on  an  explicit  characterization  of 
the  relation  between  two  Green’s  functions  and  a  tight  dimension  estimate  for  the  best 
linear  subspace  approximating  a  set  of  almost  orthogonal  vectors.  We  derive  both  lower 
bounds  and  upper  bounds  and  show  their  sharpness  for  cases  that  are  commonly  used  in 
practice. 


1.  Introduction 

Given  a  linear  differential  operator,  denoted  by  L,  the  Green’s  function,  denoted  by 
G(x,  y),  is  defined  as  the  fundamental  solution  in  a  domain  17  C  to  the  partial  differ¬ 
ential  equation 

(  L,^G(x,y)  =  (5(x-y),  x,yGl7CR" 

(1) 

I  with  boundary  condition  or  condition  at  infinity, 

where  (5(x  —  y)  is  the  Dirac  delta  function  denoting  an  impulse  source  point  at  y.  In  par¬ 
ticular,  general  solutions  of  a  partial  differential  equation  can  be  obtained  by  superposition 
of  fundamental  solutions  with  source  locations  in  17  (and/or  boundary  of  17). 

Approximate  separability  of  G(x,  y)  is  defined  as  the  following:  given  two  sets  X,Y  C 
17  C  (see  Figure  1)  and  e  >  0,  there  is  a  smallest  N'^  such  that  there  are  // (x) ,  g';  (y ) ,  /  = 
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(2) 


G'(x,y)  -'^fi{^)gi{y) 


1=1 


<  e, 


XxY 


where  ||  •  ||xxY  is  the  norm  of  some  function  space  which  G,fi,gi  belong  to.  If  X  and  Y 
are  compact  and  disjoint  domains  in  and  G(x,  y)  is  continuous  in  X  xY ,  which  is  often 
the  case  of  practical  interest,  there  exists  a  polynomial  approximation  of  G(x,  y)  in  X  x  y 
by  Weierstrass  approximation  theorem  which  is  separable.  So  there  is  a  <  oo  for  any 
e  >  0.  The  most  interesting  issue  is  how  depends  on  e,  which  manifests  the  intrinsic 
complexity  of  the  PDE  and  its  solution  within  e-approximation.  If  one  views  G(x,y)  as 
a  family  of  functions  on  X  parametrized  by  y  G  Y ,  this  is  equivalent  to  saying  that  the 
Kolmogorov  n-width  [17]  for  the  family  of  functions  G(x,y)  in  the  ||  •  ||x  normed  function 
space  is  e  when  n  =  N^.  Kolmogorov  n-width,  which  is  used  to  characterize  information 
content  in  information  theory,  is  the  best  approximation  of  a  set  5  in  a  normed  space  W 
by  a  n  dimensional  linear  subspace  defined  as 


(3) 


dn{S,W)  ;=  inf  sup  inf  \\f-g\\w, 

f£S9^L„ 


Of  course  the  role  of  x  and  y  can  be  reversed. 


Figure  1.  Green’s  function  G(x,  y)  with  dependence  on  x  G  X  and  y  G  K. 

We  introduce  the  following  relations  to  simplify  notations  in  later  derivations,  x  >  y 
means  that  there  is  a  constant  oo  >  c  >  0  such  that  x  >  cy,  x  <  y  means  that  there 
is  a  constant  oo  >  C  >  0  such  that  x  <  Cy,  and  x  ~  y  means  there  are  two  constants 
0  <  c  <  C  <  CO  such  that  cy  <  x  <  Cy.  For  our  results  for  Helmholtz  equation  (5),  all 
constants  are  independent  of  the  wave  number  k  as  k  ^  oo. 

In  this  study,  we  assume  X,  T  C  H  are  two  compact  manifolds  embedded  in  with  di¬ 
mensions  dim(X)  and  dim(y)  respectively,  i.e.,  they  may  be  compact  domains  in  R'^,  dim(X) 
dim(y)  =  d  =  1,  2,  3,  or  compact  two  dimensional  surfaces  embedded  in  R^  or  one  dimen¬ 
sional  compact  curves  embedded  in  R'^,  d  =  2,3.  Without  loss  of  generality,  we  assume 
dim(X)  >  dim(y)  =  s. 

One  can  get  a  sharper  upper  bound  for  X^  based  simply  on  regularity  of  the  Green’s 
function.  For  example,  suppose  X,  Y  are  two  disjoint  compact  domains  in  R'^  and  G(x,  y) 
is  C'™'(X  X  y),  one  can  show  that  X^  ^  using  the  following  argument.  Lay  down  a 


APPROXIMATE  SEPARABILITY  OF  GREEN’S  FUNCTION  FOR  HIGH  FREQUENCY  HELMHOLTZ  EQUATIONS 


uniform  grid  yj,j  =  1,  2, ...  J  r-^  ^  m  in  Y  with  a  grid  size  h  Vy  G  y,  G(x,  y)  can 

be  approximated  by  {m  —  l)-th  order  interpolation/extrapolation  of  G(x,  y^),  which  is  a 
linear  combination  of  G(x,  y^)  and  satisfies 

(4)  |G(-,y)  -  J]  a,G{;yj)\  <  \\D^G{;  ■)U^^x>cY)h^  <  e, 

j-yj&Bsiy) 


where  Bs{y)  is  a  ball  centered  at  y  with  a  radius  5  ^  h.  If  y  is  a  rectangular  domain 
and  G(x,  y),  one  can  even  use  We  call  G(x,y)  highly  separable  if  depends  on  e  weekly, 
W  <  0(1  loge|P),  for  some  p  >  0. 

When  a  linear  PDE,  such  as  (1),  is  discretized  and  numerically  solved,  high  separabil¬ 
ity  of  its  Green’s  function  implies  existence  of  low  rank  approximation  of  subsets  of  the 
discretized  system,  which  provides  a  matrix  compression  and  lies  at  the  heart  for  many 
efficient  algorithms.  Typically  low  rank  approximation  has  been  used  in  two  ways.  One 
way  is  to  utilize  low  rank  approximation  of  the  discretized  Green’s  function,  which  is  the 
kernel  for  boundary  integral  formulation,  for  fast  matrix  vector  multiplication  when  solv¬ 
ing  boundary  integral  equations  by  iterative  methods  [8,  11,  12,  26,  27].  Similarly,  it  has 
been  used  to  develop  fast  algorithms  for  evaluation  of  fast  oscillatory  scattering  operator 
and  Fourier  integral  operators  [6,  24,  25].  The  other  way  is  to  utilize  low  rank  property 
to  develop  fast  algorithms  to  solve  a  large  linear  system  Ax  =  b  corresponding  to  a  dis¬ 
cretization  of  a  PDE  such  as  (1).  Each  column  of  the  inverse  matrix  A~^  is  a  numerical 
approximation  of  the  Green’s  function.  Again  low  rank  approximation  for  off-diagonal 
submatrices  of  A~^,  which  is  implied  by  high  separability  of  the  Green’s  function  on  two 
disjoint  sets,  is  extensively  explored  in  many  fast  algorithms  to  solve  the  linear  system  such 
as  hierarchical  matrix  and  structured  inverse  methods  [2,  3,  7,  15,  23,  29,  31,  32].  Often 
the  low  rank  approximation  is  computed  or  learned  on  the  fly.  Fast  random  algorithms 
or  rank  revealing  type  of  methods  [14,  22]  can  be  used  to  find  the  leading  singular  values 
and  corresponding  singular  vectors  of  a  matrix.  However,  the  computation  cost  increases 
dramatically  if  the  rank  is  not  sufficiently  low.  So  both  upper  and  lower  bound  estimates 
for  approximate  separability  is  of  crucial  importance  in  these  applications. 

In  the  literature,  mostly  upper  bound  estimates  for  highly  separable  cases  were  shown 
for  Green’s  functions  or  kernel  functions  when  developing  fast  numerical  algorithms  based 
on  low  rank  approximation  as  mentioned  above.  These  estimates  are  typically  based  on 
constructive  approaches  for  Green’s  functions  or  kernel  functions  with  explicit  expression 
and  using  asymptotic  expansion.  Interesting  study  on  spatial  bandwidth  and  degree  of 
freedom  of  scattered  field  have  also  been  done  in  the  engineering  literature  [4,  5],  which 
shows  that  the  scattered  field  is  almost  band  limited  and  the  degree  of  freedom  is  close  to 
the  Nyquist  number  in  term  of  the  effective  (spatial)  bandwidth  of  the  scattered  field  and  to 
the  extension  of  the  observation  domain.  A  more  general  non-constructive  mathematical 
approach  was  developed  in  [2]  to  show  that  the  Green’s  function  for  a  coercive  elliptic 
operator  in  divergence  form  with  Loo-coefficients  is  highly  separable  (Theorem  2.8)  for  two 
disjoint  domains  X,  Y  based  on  a  key  gradient  estimate  by  Gaccioppoli  inequality.  The 
method  and  result  can  be  extended  to  Green’s  function  of  more  general  elliptic  equations 
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with  non-dominant  lower  order  terms,  such  as  convection-diffusion  equations  with  small 
convection  term  or  the  Helmholtz  equation  (5)  with  small  wave  number  k.  Their  method 
does  not  work  when  the  lower  order  term  is  dominant  which  is  the  case  for  the  Helmholtz 
equation  with  large  k.  It  becomes  a  singularly  perturbed  problem  and  the  gradient  of 
the  Green’s  function  is  unbounded  almost  everywhere  as  A:  — )■  oo  due  to  fast  oscillations. 
These  issues  are  also  reflected  in  numerical  computation  for  these  different  PDEs.  It 
is  well  known  that  there  are  many  efficient  numerical  methods  to  solve  the  discretized 
system  corresponding  to  differential  operators  that  are  elliptic  dominant,  such  as  iterative 
methods  with  various  effective  preconditioners  and  direct  inverse  methods  as  mentioned 
above.  This  is  related  to  the  intrinsic  complexity  manifested  by  the  high  separability  of 
the  corresponding  Green’s  functions.  On  the  other  hand,  it  is  well  known  that  Helmholtz 
equation  with  large  wave  number  is  very  difficult  to  solve  numerically  in  practice.  For 
example,  all  those  well  developed  iterative  methods  for  elliptic  equations  do  not  work 
effectively  for  this  case  [10]. 

Here  we  give  another  mathematics  perspective  by  showing  lower  bounds  for  the  approxi¬ 
mate  separability  of  the  Green’s  function  for  Helmholtz  equation  in  high  frequency  limit  in 
terms  of  both  e  and  k.  The  lower  bounds,  which  are  sharp  for  many  practical  setups,  show 
that  the  Green’s  function  is  not  highly  separable  as  /c  — )•  oo  and  manifests  the  intrinsic 
complexity  of  the  solution  space.  In  our  study  we  give 

•  explicit  characterization  of  the  correlation  or  angle  (in  L2  normed  space)  between 
two  Green’s  functions  of  Helmholtz  equation  (5)  in  the  high  frequency  limit, 

(||G'(-,yi)||2||G'(-,y2)||2)“Y G'(x,yi)G(x,y2)dx  <  (fcjyi  -  y2|)“",  k\yi  -  y2|  ^  oo 

Jx 

for  some  a  >  0  which  depends  on  the  dimension  of  X,  its  geometry  and  the  locations 
of  yi,y2  (see  Theorem  2.1)  based  on  generalized  stationary  phase  analysis. 

•  lower  bound  estimate  for  the  approximate  separability  for  the  Green’s  function  of 
Helmholtz  equation  in  the  high  frequency  limit 

{  a  <  f , 

Nl>\ 

{  «>§, 

and  upper  bound  estimate 

m 

as  A:  — 7-  oo  for  two  compact  manifolds  X  and  Y  with  dim(X)  >  dim(y)  =  s  and  any 
5  >  0,  where  the  constants  in  >  and  <  are  independent  of  k  for  a  fixed  small  e  (see 
Lemma  3.1  and  Theorem  3.1  -  3.4).  The  lower  bound  is  based  on  a  tight  dimension 
estimate  improved  from  that  for  a  set  of  nearly  orthogonal  random  vectors  by  N. 
Alon  [1]  and  Johnson-Lindenstraus  Lemma  [19]. 

•  explicit  estimates  and  their  sharpness  for  situations  that  are  commonly  used  in 
practice  (see  Section  4).  Our  theory  is  also  applied  to  show  precise  conditions  if 
high  separability  (or  low  rank  approximation  after  numerical  discretization)  can  or 
can  not  be  achieved  for  special  set  ups. 
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As  far  as  we  know,  lower  bound  estimate  for  approximate  separability  of  Green’s  func¬ 
tions  of  this  type  is  the  first  in  the  literature.  These  bounds  mathematically  characterize 
intrinsic  complexities  for  high  frequency  wave  phenomena.  We  hope  these  studies  and  un¬ 
derstandings  can  provide  useful  insights  for  developing  fast  numerical  algorithms  as  well. 


2.  Helmholtz  equation  and  its  Green’s  function 
Let  G{x,  y)  be  the  Green’s  function  to  the  Helmholtz  equation  in  free  space, 
(5)  AxG(x,y) -b /c^n^(x)G'(x,y)  =  (5(x  -  y), 


where  /c  >  0  is  the  wave  number,  0  <  c  <  n(x)  <  C  <  oo  is  the  index  of  refraction 
and  5{x  —  y)  denotes  a  point  source  at  y.  Suitable  far  field  radiation  condition  has  to  be 
satisfied  for  uniqueness.  The  high  frequency  limit  means  the  wave  number  /c  — )■  oo,  which 
poses  challenge  both  mathematically  and  numerically  due  to  faster  and  faster  oscillations 
in  the  solution. 

For  completeness,  we  provide  the  general  formula  for  the  free  space  Green’s  function  of 
Helmholtz  equation  (5)  for  any  space  dimension, 


(6)  G'o(x,y)  =  CrfF 


(1)/ 


ix-y|) 


d-2 


|x-y|- 


p  = 


Cd  = 


2i(27r)P’ 


x,y  G  i?'^,x  /  y. 


Hp^\r)  is  the  first  kind  Hankel  function  of  order  p  which  has  the  following  asymptotic 
behavior;  as  r  — )■  0 


2\P 


(7)  =  I  (?) 

where  r(p)  is  the  Gamma  function,  and  as  r 


oo 


p^O 

p  =  0 


(8) 


pTT 


4?  -b  0(r  2 ), 


p>  0. 


For  d  =  3,  the  Green’s  function  takes  the  simplest  form 

1  ifc|x-y| 

(9)  Go(x,y)  =  — I - r,  x/y. 

471  |x  —  y| 

For  d  =  2, 

(10)  Go(x,y)  =  -lH«(A;|x-y|)  =  x  /  y, 

and 

(11)  lim  Ho^Vr)  =  —  logr,  lim  =  \  — _j_  (Q7j,-3/2y 

r-s.0+  TT  r^oo  ^  V  vrr 
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Denote  B^{y)  and  S^(y)  to  be  a  ball  and  a  sphere  in  centered  at  y  with  radius  r 
and  p  respectively.  We  have 
(12) 

[  |G'o(x,y)p(ix  =  /  dp  f  ds  =  clujdk^^~‘^  I  r[H^^\r)fdr, 

Jo  Jsf  Jo 

d 

where  ujd  =  is  the  area  of  the  unit  sphere  in  R'^.  From  the  asymptotic  formula  (7), 
we  see  that  Gq  is  square  integrable  at  the  singular  source  for  d  =  2,3.  Also  from  the 
asymptotic  formula  (8),  we  have  ||G'o(-,  y)||2(B»(y))  ~  as  A:  — )■  oo. 

From  the  above  explicit  expressions  for  free  space  Green’s  function,  we  see  that  except 
for  the  case  d  =  3,  there  is  a  multiplication  factor  related  to  k  for  the  magnitude  of  the 
Green’s  function.  To  characterize  the  angle  or  correlation  between  two  Green’s  function 
and  study  the  separability  of  the  Green’s  function  without  the  effect  of  this  factor,  we 
define  the  normalized  Green’s  function  as 

(13)  g(X’y)  =  xeAci?'^,  ||G(-,y)||i  =  ^  |G(x,y)|2a!x, 

in  our  later  study  with  the  following  understandings: 

•  ||G(-,y)||2  is  a  smooth  function  of  y  since  fast  oscillation  due  to  rapid  change  of 
phase  function  is  removed. 

•  When  d  =  3,  all  results  for  the  normalized  Green’s  function  G{x,  y)  can  be  extended 
to  G(x,  y)  since  there  are  constants  0  <  c  <  C  <  oo  that  are  independent  of  k  such 
that  c  <  |G(x,  y)|  <  C,  Vx  G  A,  y  G  T,  once  two  compact  sets  X,Y  C  R^  are  fixed. 

•  In  this  paper  we  prove  results  for  d  =  2,3  for  practical  interest.  Since  the  Green’s 
function  is  square  integrable  at  the  source  singularity,  we  allow  overlaps  between 
two  compact  domains  X  and  Y  when  dim(A)  =  dim(y)  =  d  =  2,3.  All  results  for 
bounded  and  disjoint  X  and  Y  can  be  extended  to  d  >  3. 

Approximate  separability  of  G(x,  y)  is  defined  as  in  (2)  except  that  we  now  put  a  sub¬ 
script  k  in  to  specifically  indicate  the  dependence  on  k.  The  key  issue  is  the  dependence 
on  k  for  a  given  e  >  0  as  /c  — )■  co.  In  practice,  such  as  development  of  fast  algorithms  uti¬ 
lizing  low  rank  approximation  for  the  discretized  system,  X  and  Y  are  often  disjoint  and 
compact.  Typical  norms  used  are  either  Loo  (A  x  Y)  or  L2(A  x  Y).  In  our  study  we 
first  show  analysis  and  results  in  L2  norm,  which  fits  well  with  using  SVD  (singular  value 
decomposition)  for  low  rank  approximation  of  a  matrix,  and  then  extend  those  results  to 
Loo  norm.  Regarding  G(x,  y)  as  a  family  of  functions  in  L2(A)  parametrized  by  y  G  T 
(see  Figure  1),  the  separability  condition  (2)  in  L2(A  x  Y)  is  equivalent  to  the  existence 
of  a  linear  subspace  Sx  C  L‘2{X)  with  dimension  A|  such  that 

(14)  ^  J^\\G{x,y)  -  Ps^G{x,y)\\l^^^^dy  <  e, 

where  PsxG{x,y)  is  the  projection  of  G(x, y)  in  Sx-  This  formulation  is  the  same  with 
the  role  of  x  and  y  exchanged. 
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We  start  with  the  study  of  approximate  separability  of  the  Green’s  function,  Go(x,  y), 
for  homogeneous  medium,  i.e.,  v?{x)  =  1  in  (5),  in  free  space.  When  2-norm  is  used  as 
the  metric,  one  important  geometric  characterization  of  relation  between  two  vectors  is  the 
angle  or  correlation  between  them.  Let  X  <Z  E?  he  a.  compact  domain  and  yi,y2  ^  X  he 
two  points  with  5  =  |yi  —  y2|  <C  1.  It  is  easy  to  see  that 

1  r  e*fc(|x-yi|-|x-y2|) 

||G'o(-,yi)||2||G'o(-,y2)||2  Jx  |x-yi||x-y2| 

where  the  constant  in  <  depends  on  domain  X  and  the  distance  between  X  and  yi,y2. 
In  another  word,  two  Green’s  functions  become  more  and  more  correlated  when  the  two 
source  points  become  closer  and  closer.  This  is  true  in  general  for  Green’s  function  as 
long  as  G(x,  y)  is  Lipschitz  in  y.  Actually  for  strictly  elliptic  operator  of  the  following 
divergence  form  in  d  >  3, 


where  ajj(x)  are  bounded  measurable  functions  and  0  <  X  <  /j,  <  oo  are  two  constants, 
there  exists  a  unique  Green’s  function  [20,  13]  G{x,  y)  satisfying 

(16)  c{d,  A,  ^)|x  -  yp-*^  <  G(x,  y)  <  C{d,  A,  |u)|x  -  yp”"*, 

where  0  <  c{d,  A,  fx)  <  C{d,  A,  /u)  <  oo  are  two  constants.  Given  a  compact  domain  X  <Z 
and  two  points  yi,y2  ^  X,  define 

.  II  7^  C{d,X,i^)  r  |yi  -y2| 

p  =  mm  mm  X  -  yi  ,mm  X  -  y2  ,  K  =  ,  ,  ^ - -  1 -| - 

xsa:  xsx  c[d,X,ii)  [  p 

Then  we  have 

G'(x,y2)  ^  C{d,X,p)  'Ix-yiT  ^  C{d,X,p)  ' 

G(x,yi)  “  c{d,X,p)  l_|x-y2|J  “  c{d,X,p) 

and  vice  versa.  Given  two  disjoint  compact  domains  X,Y  C  R'^,  d  >  3,  with  p  being  the 
distance  between  the  two  domains  and  r  being  the  diameter  of  Y,  the  correlation  between 
any  two  Green’s  function  with  sources  at  yi,y2  G  T  is  bounded  by 

(17)  1  ><  G'(-,yi),G'(-,y2)  >x>  K  =  / ,  \  1  +  - 

c[d,X,p)  [  p\ 

Also  Caccioppoli  inequality  gives  a  L2  norm  bound  of  the  gradient  of  the  Green’s  function 
away  from  the  source  singularity  which  is  used  in  [2]  to  show  that  the  Green’s  function  for 
elliptic  operator  (and  Helmholtz  equation  with  small  k)  is  highly  separable.  However,  the 
picture  is  quite  different  in  the  more  challenging  regime  of  high  frequency  limit  since  the 
Green’s  function  becomes  more  and  more  oscillatory  as  /c  — )■  00. 


|x  -  y2|  +  |yi  -  y2| 
|x-y2| 


dx  —  1  <  k5 


<  Go(-,yi),Go(-,y2)  >  -1 
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2.1.  Decorrelation  of  Two  Green’s  Fnnction  in  High  Freqnency  Limit  in  Homo¬ 
geneous  Medium.  Here  we  study  the  angle  between  two  Green’s  functions  or  how  fast 
they  decorrelate  in  term  of  the  ratio  of  the  separation  distance  of  the  two  source  points 
with  respect  to  the  wave  length  due  to  fast  oscillations.  Stationary  phase  theory  will  play 
an  important  role  here.  Define 


(18) 

We  have  |cA(x)|  <  1  and 


(19) 


=  |yi  -  y2 


<^(x)  =  |yi  - 

y2|  ^(|x- 

-yil  - 

■  |x-y2|). 

V0(x)  =  |yi 

-y2r^  ( 

x-yi 

x-y2  \ 

|x-yi| 

|x-y2| ) 

r  j  (x-yi)  {: 

>^-yi)’^ 

J  (x-y2)  (x-y2)-^ 

1  ^  |x-yil 

|x-yil 

jx-y2l  lx-y2l 

|x-yi| 


|x-y2| 


|Vi;^(x)|  /  0  except  for  points  on  the  line  going  through  yi,y2  and  outside  the  interval 
between  yi,y2,  where  maximum  value  1  or  minimum  value  -1  of  (/>  is  attained  (see  Figure 
2).  Also  Il^(^(x)  is  degenerate  in  the  direction  of  yi  —  y 2-  However,  for  x  on  the  line 
and  outside  the  interval  between  yi,y2,  the  Hessian  in  the  plane  perpendicular  to  the  line, 
denoted  by  is  a  multiple  of  identity  matrix  I_\_  in  the  plane. 


(20)  Di<^(x)  =  ^ ^ 

|x-yi||x-y2| 

where  the  sign  depends  whether  maximum  or  minimum  is  attained  at  x. 

From  the  stationary  phase  result  [16,  30]  for  I{k)  =  f  where  u  G  C^{R^) 

and  (^(x)  has  isolated  stationary  points  Xm,m  =  1,2,..M,  \V(j){xm)\  =  0  and  D‘^(j){xrn) 
non-degenerate,  one  has 


pik(f)(xm) 


(21) 


11/2 


fsgn(D2</,(x™)) 


e  4 


u{Xm)\ 


where  s  >  d/2  and  C  is  an  universal  constant  independent  of  (j)  and  u.  However,  we  have  to 
modify  the  standard  stationary  phase  technique  due  to  the  following  three  complications  in 
our  case:  (1)  the  stationary  points  may  not  be  isolated,  (2)  the  integration  is  on  a  compact 
domain  X  and  the  integrand  u  is  not  C^{X),  and  (3)  there  may  be  singularities  in  the 
integrand.  Here  is  our  result. 


Theorem  2.1.  Assume  X  C  R'^,d  =  2,3  is  a  compact  domain  with  piecewise  smooth 
boundary.  Go(x,  y)  is  the  normalized  free  space  Green’s  function.  Depending  on  positions 
of  yi,y 2  relative  to  X  and  its  boundary,  there  is  some  a  >  0  such  that 


(22) 


<  Go(-,yi),Go(-,y2)  > 


^  (^|yi 


y2l)' 


min{l,  - - }  <  a  < 


d  +  1 
2 


os  k\yi  —  y2|  — )•  oo.  The  constant  in  <  depends  on  X  and  the  distances  from  yi,y2  to  X. 
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case  1 


Figure  2.  Different  positions  of  yi,y2  relative  to  X 


Proof.  We  prove  for  d  =  3  first.  Define 

k  =  k\yi  -  y2|,  (^(x)  =  |yi  -  y2rnix  -  yi|  -  |x  -  y2|), 

(23) 

—  ||Go(-,yi)||2||G'o(-,y2)||2|x-yi||x-y2|  ’ 

and  the  operator 


(24) 

we  have 

(25) 


L  = 


|V<^(x)|2^ 


dip  d 

dxj  dxj 


L'^ 


E 


i=i 


d  1  df) 

dxj  \V^{x)\‘^dxj' 


<  Go(-,yi),Go(-,y2)  > 


e*^'^Wu(x)(ix 


Denote  the  line  going  through  y i ,  y2  by  lyl  and  the  part  of  lyl  outside  the  open  interval 
between  yi,y2  by  lyl-  Depending  on  the  positions  of  yi,y2  relative  to  the  domain  X,  we 
consider  the  three  generic  cases  illustrated  in  Figure  2.  All  other  cases  can  be  deduced 
from  these  three  cases. 

Case  1.  lyl  n  X  =  0,  see  Figure  2.  Since  there  is  no  stationary  point  in  X,  i.e., 
|Vi;i)(x)|  >  c  >  0,Vx  G  X,  and  u(x)  is  smooth  in  X,  from  integration  by  part  we  have 
(26) 

eifc^(x)^(x)(ix  =  i  XY(Le*^'^W)u(x)dx 


ik 


fx  e^^^(^HL^u(x))dx  +  |V(^(x)|  2(ELi  Ui#)e*^^Wu(x)d5(x) 


Ix  e*^'^W((^'^)^«(x))dx  +  |V0(x)|  ^(ELi  uj^)e**^^WL'^u(x)o!5(x) 


+  i  Idx  I  W(x)|  Ui^)e*'='^Wu(x)o!5(x). 

Integration  by  part  can  be  continued.  However,  the  leading  term  is  the  last  term  which  is 
an  oscillatory  integral  on  the  boundary  dX.  If  the  phase  function  4>{x)  has  isolated  local 
minima  and  maxima  on  dX  and  D‘^(j){x)  is  not  degenerate  along  dX  at  those  extrema,  the 
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boundary  integral  in  the  last  term  is  <  A;  2 


by  the  stationary  phase  theory.  Hence 


(27) 


as  k  ^  00. 


If  there  is  a  piece  of  the  boundary  dX  stays  on  a  level  set  of  (/>(x),  see  Figure  3,  the  phase 
function  <^(x)  is  constant  on  that  piece,  the  boundary  integral  in  the  last  term  is  <  1  and 
hence 


(28) 


2*^^Wu(x)dx 


IX 


<k 


-1 


All  other  scenarios  are  bounded  in  between. 


as  /c  — )■  00. 


Figure  3.  A  piece  of  the  boundary  dX  stays  on  a  level  set  of  (/>(x) 

Case  2.  yi,y2  are  outside  X  but  Zyj  n  X  /  0,  see  Figure  2.  Both  <^(x)  and  u(x)  are 
smooth  in  X.  However,  all  points  on  the  line  segment  D  X  are  stationary  points  with 
the  same  phase.  Let’s  use  a  new  coordinate  system  to  evaluate  the  integral.  The  new 
orthogonal  coordinate  system  is  where  the  origin  is  at  yi  and  r-axis  is  in  the 

direction  yi  —  y2,  and  (^,  r/)  is  an  orthogonal  system  perpendicular  to  r,  see  Figure  2  (b). 
So 

(29)  f  e*^'^*'^^u(x)(ix  =  f  f 

Jx  Jri  JX{r) 

where  X{r)  denotes  the  intersection  of  X  with  the  plane  {r,^,ri)  at  a  fixed  r  and  ri  = 
L  ?’2  =  max(^^^_^)gX  r. 

For  a  fixed  r  G  [ri,  r2],  if  lyl  nX(r)  =  (r,  0, 0),  it  is  the  only  stationary  point  in  the  plane 
(r,  ^,77).  Moreover,  we  have 

4>{r,0,0)  =  1, 

(30)  Dl^ir,0,0)  = 

^l(^’®’0)lb=0  =  ||Go(-,yi)||2||Go(-,i2)||2r-(|yi-y2|+r)’ 

where  /  is  a  2  x  2  identity  matrix.  For  each  r,  one  can  use  a  partition  of  unity  for  the  domain 
X{r)  in  (^,7?)  plane:  Xi(r,  C,  h)  +  X2(r,  77)  =  l,V(r,^,77)  G  X{r).  0  <  Xi{r,^,v)  <  1  is 
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smooth  and  is  1  in  a  small  ball  centered  at  (r,  0, 0)  and  inside  X{r).  xi(^)  C;  v)  is  zero  near 

the  boundary  of  dX{r).  Denote  Ui  =  XiU,  i  =  1,  2,  then 

(31) 

Jx(r)  Jx(r)  Jx{r) 

Since  there  is  no  stationary  point  in  the  second  integral,  one  can  use  integration  by  part 
argument  as  in  case  1  to  show  that  it  is  <  k~^.  For  the  first  integral,  (r,  0,0)  is  the  only 
stationary  point.  Apply  the  standard  stationary  phase  result  and  from  (30)  we  get 


(r,  r])d^dr]  — 


l|G'o(-,yi)l|2||G'o(-,y2)||2 


<  ~k~'^ 


It  is  important  to  note  that  the  phase  in  the  leading  term  after  integration  in  (^,  rj)  over 
X(r)  is  independent  of  r,  which  means  no  fast  oscillation  when  integrating  in  r. 

For  a  fixed  r  G  [ri,r2],  if  lyl  H  X{r)  =  0,  there  is  no  stationary  point  in  X{r).  Hence 
fx(r)  rj)d^drj  will  be  less  than  the  case  when  ly\  n  X(r)  /  0.  So  we  have 


(33)  /  e*^^Wu(x)dx  =  r[ 

Jx  Jri  Jx(r) 

since  =  1  for  d  =  3. 

Case  3.  Let’s  consider  the  most  general  case  where  yi  and  (or)  are  in  the  interior 
of  X,  see  Figure  2.  The  main  contribution  still  comes  from  the  line  of  stationary  points 
ly\  n  X.  However,  singularities  of  u  at  yi  and  y2  have  to  be  taken  care  of.  Assume  that 
there  is  a  ball  with  radius  tq  <  |yi  —  y2|/4  around  each  point  yi,y2  contained  in  X. 
First,  design  a  partition  of  unity  functions,  xo(x),  Xi(x),  X2(x),  each  of  which  is  smooth 
and  non-negative  and  xo(x)  +xi(x)  +X2(x)  =  l,Vx  G  X.  Here  Xi(x),  X2(x)  are  1  in  a  ball 
centered  at  y i ,  y2  respectively  with  radius  ro/2  and  are  zeros  outside  the  ball  with  radius 
ro.  xo(x)  =  1  -  xi(x)  -  X2(x).  Denote 

u(x)  =  u(x)  [xo (x)  +  Xi  (x)  +  X2 (x)]  =  Uq (x)  -k  Ul  (x)  +  U2 (x)  . 

We  break  the  integral  in  (25)  into  three  parts: 

(34)  f  e*^'^*'^^u(x)dx  =  f  e*^^*-’’‘^(uo(x)  +  ui(x)  -|-  U2(x))dx  =  I-I-H+HI. 

Jx  Jx 

The  first  term  can  be  reduced  to  Case  2.  Now  let’s  look  at  the  second  term  in  (34).  We 
change  the  integration  to  a  spherical  coordinate  (r,  9,  V’)  centered  at  yi  with  6  G  [0,  27r] 
being  the  azimuthal  angle,  V’  £  [0;  tt]  being  the  polar  angle  and  yi  —  y2  being  the  polar 
axis.  Then 

(35)  [  e*^^Wui(x)dx=  r  [  e*^'^Wui(x)dsdr 

Jx  Jo  JdB{yi,r) 


12 


BJORN  ENGQUIST  AND  HONGKAI  ZHAO 


It  can  be  seen  from  (19)  that  V(/>(x;  yi,  y2)  is  never  aligned  with  the  normal  at  x  of  the 
sphere  centered  at  y2  except  at  the  intersections  of  lyl  with  the  sphere.  So  on  any  sphere 
dB{yi,r)  there  are  exactly  two  stationary  points  at  'ip  =  0  and  ip  =  tt  which  are  non¬ 
degenerate.  At  the  two  stationary  points  we  have 

(p{r,9,0)  =  1 


(36) 


n  n\  ^  _ Xlir,0,0) _ 

^{r,e,7r)  =  |yi  -  y2\~^{‘^r  -  |yi  -  y2|) 


Dl4>{r,  e,Tr) 


|yi-y2|-2r  j- 

»’{|yi-y2|-0 


tti(r,  0,71) 


_ Xlir,d,TT) _ 

||Go(-,yi)||2||Go(-,y2)||2r-(|yi-y2|-r) 


where  _L  denotes  the  tangent  plane  of  the  sphere.  Note  that  modulo  a  scaling  factor  r 
ui  and  its  derivatives,  and  D'^cp  as  functions  on  dB{yi,r)  are  all  smooth  and  uniformly 
bounded,  i.e.,  \D^ui\  =  0{r~^)  and  ||Z)^(^||  =  0{r~^).  After  scaling  dB{yi,r)  to  the  unit 
sphere  and  apply  the  stationary  phase  result  (21)  to  the  two  stationary  phase  points,  one 
gets 


(37) 


l|Go(-,yi)||2||Go(-,y2)||2 


-e 


ik 


Xi(r,6»,0)  + 


e^fc|yi-y2l  ^(2^-|yi-y2l);^]^(r,0,7r)- 
|yi-y2|-2r 


2  ^ 


The  righthand  side  in  the  above  expression  comes  from  an  estimate  of  the  righthand  side 
term  of  the  stationary  phase  formula  (21)  and  the  constant  in  <  is  uniformly  bounded 
when  r  — 7-  0.  The  first  term  in  the  square  bracket  is  the  leading  term  from  the  stationary 
phase  at  =  0  on  the  sphere  dB{yi,r)  and  the  phase  is  constant  in  r.  The  second  term 
in  the  bracket  is  the  leading  term  from  the  stationary  phase  at  V’  =  ^  on  the  sphere 
dB{yi,r).  However,  it  has  a  phase  dependent  on  r  which  results  in  a  higher  order  term 
after  integration  in  r.  Since  all  terms  are  integrable  as  r  — )•  0,  we  have 


(38) 


rro  r  -  - 

/  e*^‘^Wui(x)dx 

= 

/  /  e*^'^Wui(x)dsdr 

Jx 

Jo  J  dB{yi,r) 

The  third  term  in  (34)  can  be  shown  in  the  same  way. 

From  the  above  analysis  we  see  that  the  main  contribution  for  the  integral  of  <G'o(')  yi),  Go(')  y2)>A' 
may  come  from  the  stationary  line  and/or  the  boundary  integral  on  after  integration 
by  part.  All  other  cases  can  be  reduced  to  the  above  three  cases. 
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In  2D,  the  Green’s  function  in  free  space  is  of  the  form  (10)  with  the  asymptotic  formulas 
(11).  For  Case  1  and  2,  the  asymptotic  formula  (11)  for  r  — )•  oo  can  be  used  as  A:  — )■  oo. 

Since  the  phase  function  in  the  exponential  for  <  G'o(-,  yi),  G'o(-,  y2)  >  is  also  of  the  form 
A;(|x  — yi|  —  |x  — y2|),  same  arguments  used  above  can  be  applied.  So  we  have  the  following 
analogous  results  in  2D: 

Case  1.  Since  the  boundary  dX  is  a  one  dimensional  curve,  there  is  some  a,  1  <  a  < 

=  I ,  such  that 

(39)  <  Go(-,yi),G'o(-,y2)  >  <  (A:|yi  -  y2|)“",  as /c|yi  -  y2|  ^  oo 

Case  2.  The  leading  contribution  is  due  to  the  line  of  stationary  phase  ly\  except  that  the 
dimension  orthogonal  to  the  line  is  ID,  we  have  |  and 

(40)  |<  G'o(-,yi),Go(-,y2)  >1  <  (A:|yi  -  y2|)"T  as  A:|yi  -  y2|  oo 


Case  3.  The  singularity  at  the  source  is  also  integrable  hence 

(41)  |<  G'o(-,yi),Go(-,y2)  >1  <  (A:|yi  -  y2|)"T  as  A:|yi  -  y2|  oo 


□ 


Here  we  give  a  few  remarks  related  to  the  theorem  above. 

Remark  2.1.  For  d  =  3,  the  same  estimate  also  holds  for  two  unnormalized  Green’s  func¬ 
tions,  i.e.,  \<  Go(-,yi),Go(-,y2)  >|  ~  <  Go(-, yi),  G'o(-, y2)  >  ,  since  0  <  c  <  ||G'o(-, y)||L2(A)  < 
C  <  CO  as  k  ^  oo  for  two  constants  c  and  C  that  only  depend  on  X .  However,  this  is  not 
true  ford  =  2.  // yi,  y2  ^  X,  |<  Go(-,  yi),  Go(-,  y2)  >|  ~  <  G'o(-,  yi),  Go(-,  y2)  >  as 

k  ^  oo  due  to  the  asymptotic  formula  (11). 

Remark  2.2.  The  estimate  in  Theorem  2.1  characterizes  the  correlation  or  angle  between 
two  normalized  Green’s  function  in  term  of  the  ratio  of  the  separation  distance  between  the 
two  sources  and  the  wavelength.  One  can  also  incorporate  another  scaling  factor  when  the 
distance  from  the  two  points  y*,  i  =  1,  2  to  X  is  large  compared  to  yi  —  y2| .  Geometrically, 
this  means  that  the  difference  between  two  distance  functions,  |x  —  yi|  —  |x  —  y2|,  changes 
slowly  with  respect  to  x  £  X.  Hence  fast  oscillation  due  to  rapid  change  of  the  phase 
function,  i/c(|x  —  yi|  —  |x  —  y2|)  is  discounted  and  the  decorrelation  rate  of  two  Green’s 
function  is  reduced. 

Assume  the  size  of  X  is  0(1)  (otherwise  one  can  first  scale  x  by  the  size  of  X  for  the 
Helmholtz  equation  (5) )  and  ~  dlsdyf  x)  ~  ^  1;  which  falls  into  either  Gase  1 

or  Gase  2  in  Theorem  2.1.  From  (19),  we  see  that  Vf  is  scaled  by  p  and  are 

scaled  by  p,  p'^  respectively  when  they  are  not  degenerate.  The  scaling  for  u(x)  is  canceled 
due  to  the  normalization  according  to  the  definition  (23).  When  applying  the  stationary 


14 


BJORN  ENGQUIST  AND  HONGKAI  ZHAO 


phase  result  (21)  at  a  point  of  stationary  phase,  one  can  see  that  k  is  rescaled  to  kp  if 
kp  ^  oo.  For  Case  1,  the  main  contribution  for  <  Go(',  yi),  Go(-,  y2)  >  comes  from  the 

last  term  in  (26).  There  is  a  scaling  factor  of  p~^  from  |Vi^|“^  due  to  the  integration  by  part 
and  there  are  isolated  stationary  points  on  dX  in  general.  So  overall  k  is  rescaled  to  kp  for 
Case  1.  For  Case  2  in  Theorem  2.1,  the  main  contribution  for  <  Go(-,  yi),  Go(-,  y2)  > 

comes  from  the  line  of  stationary  phase  in  general,  where  Df  and  are  degenerate. 
Along  the  line  of  stationary  phase,  and  det[D'j_f]  is  scaled  by  p^  and  respec¬ 

tively  from  (20).  Applying  the  stationary  phase  result  (21)  in  the  plane  perpendicular  to 
the  line  of  stationary  phase,  k  is  rescaled  to  kp^.  From  both  cases  we  see  that  k  is  at  least 
rescaled  to  kp  in  the  decorrelation  estimate  for  two  Creen’s  function. 

Remark  2.3.  One  can  generalize  the  arguments  in  Theorem  2.1  to  more  general  situations 
where  X  is  a  compact  manifold  embedded  in  with  dim(X)  =  s  <  d,  such  as  a  surface 
(s  =  2)  in  or  a  curve  (s  =  1)  in  d  =  2,3.  For  example,  if  X  is  a  compact  manifold 
without  boundary,  e.g.,  a  closed  surface  or  curve,  and  two  points  yi,y2  ^  X,  there  is  some 
a,  0  <  a  <  ^ 

(42)  <  Go(-,yi),Go(-,y2)  >  <  (A;|yi -y2|)"",  as  A:|yi  -  y2|  ^  oo. 

The  two  extreme  cases  are:  (1)  a  =  0  happens  when  there  is  a  piece  of  X  stays  on  a  level 
set  of  the  phase  function  (/>(x)  =  |yi  —  y2|“^(|x  —  yi|  —  |x  —  y2|);  (2)  a  =  ^  happens  if  the 
phase  function  0(x),  which  has  stationary  phase  points  on  a  compact  manifold,  has  isolated 
stationary  phase  points  on  X  and  is  not  degenerate  along  X  at  those  points.  The 

later  case  is  more  generic. 

If  X  is  a  compact  manifold  with  boundary,  there  is  some  a,  0  <  a  <  such  that 

(43)  <  Go(-,yi),Go(-,y2)  >  <  (A;|yi -y2|)"",  as  A:|yi  -  y2|  ^  oo. 

The  two  extreme  cases  are:  (1)  a  =  0  happens  when  there  is  a  piece  of  X  stays  on  a  level 
set  of  the  phase  function  (2)  a  =  happens  if  the  phase  function  0(x)  has  no 

stationary  phase  in  X  and  has  isolated  stationary  phase  points  on  dX  and  Z1^0(x)  is  not 
degenerate  along  dX  at  those  points.  If  there  are  isolated  stationary  phase  points  in  the 
interior  of  X,  a  =  ^. 

Remark  2.4.  According  to  the  Hessian  estimate  (20),  there  are  two  axisymmetric  k  de¬ 
pendent  domains  around  the  stationary  line  ly\  on  each  side  of  yi  and  y2,  denoted  by  Ri 
and  i?2  respectively  (see  Figure  4),  in  which  the  phase  function  kf)  does  not  change  rapidly. 
For  example,  let’s  look  at  a  point  x  G  lyl  on  the  side  of  y2  and  denote  r  =  ±|x  —  y2|. 
Again  we  use  the  coordinate  system  x  =  (r,  ^,ry)  as  shown  in  Figure  4-  Since  (/)(r,  0, 0)  = 
1,  |V(^(r,  0, 0)1  =  0,r  >  0,  for  a  point  {r,f,ri)  with  r  >  0,  ,  from 
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(20)  we  have 
(44) 


fe|V(^(r,^,7?)|  < 


fe|yi  -y2|\/C^  + 

r{r  +  |yi  -  y2|) 


<  1. 


Hence  <  Go{-,yi),Go{-,y2)  >  is  not  an  oscillatory  integral  in  Ri  or  R2  and  the  two  Green’s 
functions  do  not  decorrelate  fast  in  a  subdomain  contained  in  these  two  domains.  We  will 
also  provide  special  k  dependent  setups  of  domains  X,Y  such  that  G(x,  yi)  and  G(x,  y2) 
do  not  decorrelate  fast  in  X  for  two  sources  yi,y2  £  Y  in  Section  4-2- 


Figure  4.  A  domain  where  two  Green’s  functions  do  not  decorrelate  fast. 


Remark  2.5.  The  correlation  between  two  Green’s  function  can  also  be  used  for  study 
imaging  resolution  using  waves.  Suppose  X  is  a  compact  planar  region  in  R^  where  the  wave 
field  is  measured.  yi,y2  G  R^  are  two  point  sources  or  scatterers.  If  the  line  connecting  yi 

3 

and  y2  is  parallel  to  X,  we  have  \<  Go(-, yi),  Go(-,  y2)  >|  ^  (^lyi  ~  y2\)~^  in  general  as 
^|yi  “  y2|  — >  00  since  there  is  no  stationary  phase.  While  if  the  line  connecting  yi  and  y2 
intersects  X  perpendicularly,  |<  G'o(-, yi), Go(-, y2)  >|  <  (fe|yi-y2|)“2  asA:|yi-y2|  -^00 
since  the  intersection  point  is  a  stationary  point.  Hence  it  implies  imaging  resolution  in 
the  range  direction  is  poorer  than  that  in  the  plane  parallel  to  X . 

2.2.  Decorrelation  of  Two  Green’s  Function  in  Inhomogeneous  Medium  in  High 
Frequency  Limit.  The  situation  is  similar  for  inhomogeneous  medium  in  the  high  fre¬ 
quency  limit  when  the  Green’s  function  has  a  valid  geometric  optics  approximation.  The 
assumption  for  geometric  optics  ansatz  is  that  the  solution  to  the  Helmholtz  equation  has 
the  following  expansion: 

00 

(45)  u(x)  =  a^(x)(ife)-™, 

m=0 

where  <^(x)  is  the  phase  function  and  am(x)  are  the  amplitude  functions  which  satisfy 

|V(/)(x)|  =  n(x),  ■  Vao(x)  -|-  A(/)(x)ao(x)  =  0 

2V(/)(x)  •  Vam(x)  A(/)(x)am(x)  Aam-i(x)  =  0,  m  =  1,  2, . . . 


(46) 
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For  the  Green’s  function  G(x,  y)  defined  in  (5),  we  have  the  following  condition  at  the 
source  point  y: 

(47)  lim  I  -R(y)]  =  0,  lim  ao(x,  y)47r|x  -  y|  =  1. 

x-5>y  \  |x  —  y|  J  x-?>y 

The  above  geometric  optics  ansatz  can  be  formulated  as  a  Hamiltonian  system,  also 
called  as  Lagrangian  formulation  or  ray  tracing,  which  gives  explicit  ordinary  differen¬ 
tial  equations  (ODE)  for  bicharacteristics  (x(t),p(t))  in  phase  space  with  Hamiltonian 
H (x,  p)  =  IpI  —  n(x)  and  p  =  Vcj), 


dx(f,xo,po) 

dt 


Vp77(x,p) 


P 

n(x)  ’ 


x(0)  =  xo, 


(48) 


dp(f,xo,po) 

dt 


-Vxi7(x,p)  =  Vn(x),  p(0)  =  po  =  V(/)(xo) 


#(x(t,xo,po)) 

dt 


n(x), 


</>(0)  =  </>(xo) 


The  projection  of  the  bicharacteristics  in  the  physical  space,  i.e.,  x(t,  yo,po),  are  called 
rays.  If  there  is  no  caustics,  i.e.,  two  rays  do  not  intersect  in  the  physical  space,  each  ray 
is  a  geodesic  in  the  physical  space  with  the  slowness  n{x)  =  as  the  metric,  where  c(x) 
is  the  wave  speed.  |0(x(t2,  xq,  po))  —  (/>(x(ti,xo,  po))|  is  the  shortest  travel  time  between 
points  x(ti,xo,po)  and  x(t2,xo,po).  Moreover,  the  amplitude  along  each  ray  is  given  by 


(49) 


ao(x(t,xo))  =  ao(xo) 


n(xo) 

n(x(t,xo)) 


9x(t,xo) 

9X0 


where 


9x(t,xo) 

0X0 


is  the  determinant  of  the  Jacobian 


0x(l,Xo) 

0X() 


meaning  the  geometric  spreading 


of  rays.  Before  caustics  are  formed,  the  determinant  is  always  positive  and  bounded.  Once 
ao(x)  is  known,  ai(x),  02  (x), . . .  can  be  solved  consecutively  from  (46). 

In  particular,  for  the  geometric  optics  ansatz  for  the  Green’s  function  with  a  point 
source  at  yo,  rays  x{t,yo,9)  are  emanating  from  the  source  yo  in  all  directions  and  can 
be  parametrized  by  the  initial  directions,  i.e.,  the  take-off  angles  0  G  5"^“^  on  a  unit 
sphere.  The  ODEs  for  the  rays  (48)  have  initial  conditions  x(O,yo,0)  =  yO)P(0,  yo,^)  = 
n{yo)0,4’{O,yo,6)  =  0,  with  ao(O,0)  evenly  distributed  in  all  6.  If  there  is  no  caus¬ 
tics,  every  point  x  has  a  unique  ray  passing  through  it,  i.e.,  Vx,  3!t(x),  0(x)  such  that 
x(t(x),  yo,  0{x))  =  X.  Moreover,  the  ray  connecting  yo  and  x  is  the  geodesic,  or  the  short¬ 
est  travel  time,  between  the  two  points  with  n(x)  being  the  slowness  of  the  medium.  Under 
no  caustics  assumptions,  the  Greens  functions  with  sources  at  yi,y2  can  be  approximated 
by  the  following  geometric  optics  ansatz 


(50) 


G(x,y,)-e*"'^(-’y9yl^.(x)  <  j  =  1,2, 


APPROXIMATE  SEPARABILITY  OF  GREEN’S  FUNCTION  FOR  HIGH  FREQUENCY  HELMHOLTZ  EQUATIONS 


where  Aj(x)  =  J2m=o^rn{'^,yj){ik)  ™-,j  =  1,2.  Hence 


(51)  <  G(-,yi),G(-,y2)  >  -  / 

Jx 

where 

k  =  k4>{yi,y2),  (^(x)  =  (/>"^(yi,y2)(</>(x,yi) -(/)(x,y2)),  u(x)  = 


Hi(x)yl2(x) 


||G(-,yi)||2||G(.,y2)||2 


Denote  Fy^  to  be  the  unique  ray  that  passes  through  yi  and  y2  as  illustrated  in  Figure  5 
(a).  If  n(x)  is  smooth  and  0  <  c  <  n(x)  <  C  <  oo,  one  has  (see  Figure  5(b)) 


c|yi  -  y2|  <  </>(yi,  y2)  <  G|yi  -  y2| 


(x,yi)  -  (/>(x,y2)|  <  (/)(y2,yi) 


So  the  phase  function  (/>(x)  attains  the  global  maximum  or  minimum  ±1  on  the  part  of 
the  ray  Fy^  which  is  outside  the  interval  between  yi  and  y2,  denoted  by  Fy^.  Moreover, 
Vx<^(x,  yi)  —  Vx</>(x,  y2)  /  0  for  any  x  that  is  not  on  Ty\  because  the  two  different  and 
unique  geodesics  connecting  x,  yi  and  x,  y2  respectively  can  not  be  tangent  at  x  .  So  Fyj 
is  a  stationary  curve  in  the  inhomogeneous  case  and  plays  the  same  role  as  the  straight 
line  lyl  in  the  homogeneous  case.  Depending  on  whether  the  ray  passes  through  X  or  not, 
we  get  the  same  results  as  in  Theorem  2.1  because  A:(/>(yi,y2)  ~  k\yi  —  y2|.  This  is  true 
when  yi  and/or  y2  are  in  X  as  well  since  the  amplitude  uq  satisfying  (47)  has  exactly  the 
same  singularity  as  the  homogeneous  Green’s  function. 

The  main  complication  for  geometric  optics  ansatz  in  heterogeneous  medium  is  when 
rays  cross  each  other,  i.e.,  when  caustics  are  formed.  Although  bicharacteristics  in  phase 
space  are  still  well  defined,  the  amplitude  formula  (49)  breaks  down.  So  in  general,  the 
above  arguments  can  not  be  carried  over  to  a  general  inhomogeneous  medium.  However,  in 
the  case  that  there  are  finite  number  of  rays  starting  from  yi  going  through  y2  and  there  is 
a  partition  of  unity  for  the  takeoff  angle  6  on  such  that  there  is  a  small  cone  around 
each  ray  where  there  is  no  caustics,  then  one  can  apply  the  above  argument  to  each  cone 
and  get  the  same  results. 


Figure  5.  Rays  in  inhomogeneous  medium 
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3.  Approximate  Separability  Estimate  for  the  Green’s  function  of 
Helmholtz  equation  in  High  Frequency  Limit 

In  this  section  we  present  general  estimates  for  the  approximate  separability  of  Green’s 
function  of  Helmholtz  equation  in  the  high  frequency  limit.  In  the  next  section,  we  will 
apply  these  results  to  get  explicit  estimates  for  special  setups  that  are  of  interest  in  practice. 
These  estimates  imply  rank  estimates  for  discretized  operators  which  are  important  for 
developing  fast  algorithms  for  solving  high  frequency  Helmholtz  equation  and  its  boundary 
integral  equation  counterpart. 


3.1.  Approximating  a  Set  of  Almost  Orthogonal  Vectors  by  a  Linear  Subspace. 

First  we  present  some  background  and  introduce  definitions  for  the  approximation  of  a  set 
of  vectors  using  a  linear  subspace.  It  will  be  extended  later  to  the  approximation  of  Green’s 
function  in  the  infinite  dimensional  function  space.  Let  G  R'^,  m  =  1,2, . . . ,  N  he  a  set 
of  vectors.  Define  matrix  V  =  [vi,  V2, . . . ,  vw]  and  matrix  A  =  [amn]NxN  =  V'^V.  Let 
Ai  >  A2  >  . . .  >  Aat  >  0  be  the  eigenvalues  of  A,  then  tr{A)  =  J2m=i  =  Ylm  Il’'^m|l2- 
\/A7  >  \/A2  >  • . .  >  V^N  ^  0  are  also  called  singular  values  for  V.  The  best  linear 
subspace  Si  of  all  linear  subspace  Si  of  dimension  I  that  approximates  the  set  of  vectors 
{vm}m=i  ™  least  square  sense  is  the  space  spanned  by  the  first  I  left  singular  vectors  of  V 
and  satisfies 

N  N 

(52)  y'  l|Vm  -  Vmlli  =  min  y]  ||Vm  -  PSi^mWl 

S]  ,dini(5')=/ 

m  ‘  ^  m 


N 


E 


where  Ts'jV  denotes  projection  of  v  in  Si.  In  another  word, 


(53) 


A;  = 


eeR‘*,| 


max 

e||2  =  l,e±S;_l 


N 

E 

j 


Vo 


is  the  maximum  reduction  of  approximation  error  in  term  of  least  square  for  the  set  of 
vectors  {vm}m=i  when  adding  one  more  dimension  to  the  previous  optimal  1  —  1  dimensional 
linear  subspace.  Here  we  introduce  two  definitions  for  approximate  rank  estimate  for  a 
symmetric  non- negative  matrix  A. 

Definition  3.1.  Given  e  >  0,  =  maxi<m<Ar  RJ)  s.t.  vXn  >  e,  i.e.,N^  denotes  the 

largest  m  such  that  \/Xm  >  e. 

Definition  3.2.  Given  1  >  e  >  0,  =  minM,  s.t.  Ylm=M+i  Am  ^  Ylm=i  Am- 

If  A  =  V’^V,  definition  3.2  implies  that  if  a  linear  subspace  S'"  satisfies 


E 


II 

'N 

m=l 


-  -Pg^Vmlli 

llv  l|2 
II  ^'m\\2 


<  e 


2 


(54) 
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then  dim(5'^)  >  Nf.  Assume  0  <  c  <  ||vm||2  <  C  <  oo,  m  =  1,  2, . . . ,  A",  we  can  conclude 
that  if  a  linear  subspace  satisfies 


Em=l  l|Vm  - 

N 


(55) 


<  ce 


E 


N 


-  Ps^^n 


Z-/m=l 


< 


then  dim(S"^)  >  iV*^.  In  another  word,  the  least  dimension  of  a  linear  subspace  that  can  have 
an  ce-r.m.s.  (root  mean  square)  approximation  of  a  set  of  vectors  v^,  m  =  1,2, . . . ,  N  is 
at  least  N^.  The  root  mean  square  approximation  will  lead  to  L2  approximate  separability 
estimate  for  Green’s  function  in  the  continuous  case. 

In  the  previous  section,  we  proved  the  rate  of  decorrelation  of  two  Green’s  function:  |  < 
G'(-,yi),G'(-,yi)  >  |  <  (/c|yi  -y2|)“"  for  some  a  >  0  as  k\yi  -y2\  00.  Geometrically  it 

means  that  two  Green’s  functions  with  sources  separated  a  little  more  than  one  wavelength 
become  almost  orthogonal  as  /c  — )■  00.  Intuitively,  for  two  domains  X,Y  C  if  one  views 
G(x,  y)  as  a  family  of  functions  in  L2{X)  parametrized  hy  y  gY  and  lays  downs  a  uniform 
grid  yj  G  Y  with  grid  size  h  =  k~^  for  any  0  <  /3  <  1,  G{'k,  yj)  is  a  set  of  almost  orthogonal 
vectors  in  L2{X).  A  natural  question  is  the  least  number  of  dimensions  of  a  linear  space 
that  can  contain  a  set  of  almost  orthogonal  vectors.  This  question  has  been  studied  in  [1]  by 
rank  estimate  for  small  off-diagonal  perturbation  of  identity  matrices,  which  is  equivalent  to 
the  same  question  for  a  set  of  almost  orthogonal  unit  vectors.  In  particular,  the  asymptotic 
estimate  is  optimal  and  is  used  to  show  the  sharpness  of  Johnson-Lindenstraus  Lemma  [19]. 
However,  this  result  can  not  address  our  problem  adequately  for  the  following  two  reasons. 
First,  the  assumption  in  [1]  on  almost  orthogonality  for  a  set  vectors  is  only  pairwise. 
In  our  problems,  the  set  of  vectors  are  Green’s  functions  for  a  PDF  which  has  spatial 
structure,  i.e.,  the  angle  between  two  Green’s  functions  depends  on  separation  distance  of 
the  two  sources.  The  spatial  structure  has  to  be  taken  into  account  to  get  sharp  estimates. 
Second,  approximate  separability  means  that  we  need  to  estimate  the  least  dimension  of 
a  linear  subspace  that  can  approximate  a  set  of  vectors  to  a  certain  tolerance  instead  of 
containing  the  whole  set  of  vectors.  Here  we  adopt  the  approach  from  [1]  and  develop 
more  careful  estimates  in  Lemma  3.1  for  a  set  of  Green’s  functions  by  taking  into  account 
both  the  spatial  structure  and  approximation  tolerance.  The  lemma  is  then  used  to  prove 
lower  bound  estimates  for  approximate  separability  of  the  Green’s  function  for  Helmholtz 
equation  in  the  high  frequency  limit. 

For  more  generality,  we  assume  X,Y  are  two  compact  manifolds  embedded  in  R^,  i.e., 
they  may  be  compact  domains  in  R^,  or  compact  surfaces  embedded  in  R^  or  compact 
curves  embedded  in  R'^,  d  =  2  or  3.  Without  loss  of  generality,  we  assume  d  >  dim(X)  > 
dim(y)  =  s.  For  a  smooth  manifold  Y  with  dim(y)  =  s  =  1,2,3,  it  contains  a  local 
patch  of  size  0(1)  that  is  diffeomorphism  to  a  one  dimensional  line  of  unit  length,  a  two 
dimensional  unit  square  and  a  three  dimensional  unit  cube  respectively,  which  we  will  only 
consider  in  the  following  analysis.  Let  G(x,  y)  be  the  Green’s  function  of  the  Helmholtz 
equation  (5).  The  following  Lemma  3.1  shows  the  dimension  estimate  for  a  linear  subspace 
in  L2{X)  which  approximate  of  a  set  of  Green’s  function  G(x,  ym),ym  G  Y  with  e-r.m.s 
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Lemma  3.1.  Let  X,Y  be  two  compact  manifolds  embedded  in  R'^,d  =  2,3  and  d  > 
dim(X)  >  dim(y)  =  s.  If  for  any  two  points  yi,y2  G  Y, 

I  <  G(-,yi),G(-,y2)  >  I  <  {k\yi  -y2|)"“  as  /c|yi  -y2|  ^  oo 
for  some  a  >  0,  then  there  are  points  ym  &  Y,m  =  1,2, ,  iV|  ~  ,  for  any  0  <  (5  <  1 

TV’S 

and  arbitrary  close  to  0,  such  that  for  the  set  of  Green’s  functions  {G{x,ym)}m=i  I^2{X) 
and  matrix  A  =<  G{-,ym),  G{-,yn)  > 

f  (l-e2)2fc2-,  a<|, 


Nl  > 

—k 


{l-e‘^fk^-\  «>§, 


Nl  < 


-4:j^2{s—a-5) 


..-Aus-S 


a  <  2, 


a  >  2) 


as  k  ^  oo,  where  the  constants  in  <  and  >  only  depend  on  X,  Y  and  n(x). 

Proof.  We  prove  the  statement  for  X,Y  <Z  B?  and  dim(y)  =  s  =  1,2,3  respectively.  The 
case  for  X,Y  <Z  P?  can  be  proved  in  exactly  the  same  way. 

Case  1;  s  =  1.  T  is  a  line  of  unit  length  in  P? .  Put  down  a  uniform  grid  ym  in  Y  with 
grid  size  h  =  k^~^ ,  0  <  d  <  1,  such  that  \ym  —  yn|  ~  |m  —  nlh,  m,n  =  1,2, ...  ,n^  =  k^~^ . 
See  Figure  6(a).  Dehne  the  matrix 


where 


A  —  {Umn )«'*  X  ’  Rmn  —  ^  Cr(.' tY m)  7  G(^- ,y ji)  P  7 

k  k 


G>mm  —  1;  |o!"m,n|  ^  '^1  k  ,  ?Tl,  71  —  1,2,...,  Tlf.. 


Let  Ai  >  A2  >  . . .  >  A^/i  >  0  be  the  eigenvalues  of  A.  Then  Since  = 

''k 

m.&y.2<m<n>l  ^  e  and  =  minM,  S.t.  YTm=M+l  ^rn  <  Y2m=l  ^rn  = 

we  have 


Am  >  (1  -  e^)  ^  Am  =  (1  -  e^)n'l 


E  >  E  ^ 


(61) 
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At  the  same  time,  for  a  fixed  a  >  0  and  take  0  <  5  <  1  arbitrary  close  to  0, 


E 


n 


h 

k 


m=l 


E 


n 


h 

k 


n=l 


a 


2 

mn 


(62) 


2  I  o  ^ 

^mm  ^  Y2m=n-\-l  ® 


'  _|_  ^2(l-(5-a)  < 

<<  In/c  < 

_  +  ^l-5-2a5  <  ^1-5^ 


2 

m,m—n 


a  <  ^, 
a  =  i 

Cl  2 ) 

a  >  5, 


2a  <  1  -  5  <  1 
0  <  (5  <  1 
0  <  (5  <  1 


Hence  for  a  fixed  a  >  0  and  any  0  <  5  <  1  arbitrary  close  to  0,  combining  (60)  and  (62) 
we  have 

r  (1  -  a  <  I 

(63)  iV^  >  <^ 

[  (l-e2)2^i-5^  a>i 

and  combining  (61)  and  (62)  we  have 

(  Ol  <\ 

(64)  Nl<{ 

i  e  a  >  5 

Case  2;  s  =  2.  H  is  a  unit  square  in  R^.  Again  put  down  a  uniform  grid  ym-,'ni  = 
1,  2, . . . ,  =  ^2(1-5) 

in  Y  with  grid  size  h  =  ^,0<(5<1  (see  Figure  6(b)),  and  define 

matrix  A  as  in  (58).  Let  Ai  >  A2  >  •  •  •  >  >  0  be  its  eigenvalues,  then  Ylm=i  =  nt 

^k 

and  we  have  (60),  (61).  At  the  same  time 


/I'k  '^k  '^k 

(65)  51  =  tr{A^A)  =  51  51 

ra=l  m=l  n=l 

Let’s  look  at  the  sum  of  each  row.  Assume  is  the  center  of  the  square.  We  divide  all 
other  points  into  groups  of  1st  square  neighbors,  2nd  square  neighbors,  . . .,  j-th  square 
neighbors,  denoted  by  5j,  j  =  1,  2, . . . ,  J  ~  h~^  =  kX~^ .  See  Figure  6(b).  Sj  contains  those 
4(2j  +  1)  —  4  =  8j  grid  points  that  are  on  the  4  sides  of  the  square  centered  at  with 
each  side  of  length  2jh.  Then  we  have  jh  <  jy^^  —  ym|  <  V2jh,ynj  G  Sj,  and  hence 
0'm,nj  ^  {kjh)~°^  =  j~°‘k~°‘^ .  For  a  fixed  a  >  0  and  take  0  <  <5  <  1  arbitrary  close  to  0,  we 
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have 


(66) 


En=l  «mn  =  1  +  E/=1  ^  1  +  E/=1  J 

'  1  +  a<l,  a<l-(5<l 

<  <  1  +  ln/c<l,  a  =  l,  0<(5<1 

l  +  <1^  a>l,  0<(5<1 


-2a5 


Actually  for  any  grid  point  y^,  each  j-th  square  neighbors  of  has  at  least  8j'/4  =  2j 
points  for  j  =  1,  2, . . . ,  J  ~  k^~^,  e.g.,  if  y^  is  a  corner  point  of  the  square  domain.  Hence 
the  above  asymptotic  formula  is  still  true.  For  a  fixed  a  >  0  and  any  0  <  (5  <  1  arbitrary 
close  to  0,  we  have 


(67) 


m=l 


X]  Z]  ~ 

m=l  n=l 


fe2(2(i-<5)-«),  a<l 
a  >  1 


Combining  (67)  with  (60),  we  have 


(68)  iV|  >  I 

Combining  (67)  with  (61),  we  have 


(l_g2)2^2a^ 

(l_e2)2p(i-5)^ 


(69) 


Nl  < 


g-4p(2(l-5)-a)^ 

,-4^2(l-J)^ 


a  <  1 
a  >  1. 

a  <1 
a  >  1 


yi 

yi 


Y 


(a) 


Figure  6.  Green’s  functions  with  sources  on  a  uniform  grid. 
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Case  3:  s  =  3.  y  is  a  unit  cube  in  R^.  Put  down  a  uniform  grid  ym,  m  =  1,2, . . .  ,n^  = 
p(l-5) 

in  Y  with  grid  size  h  =  k^  ^,0  <6  <1  and  define  matrix  A  as  in  (58).  Let 

1 

Ai  >  A2  >  • . .  >  X^h  >  0  be  its  eigenvalues,  then  X]rn=i  Am  =  and  we  have  (60), 

k 

(61).  Similar  to  2D  case,  assume  y™  is  the  center  of  the  cube.  We  divide  all  other  points 
into  groups  of  1st  cube  neighbors,  2nd  cube  neighbors,  . . .,  j-th  cube  neighbors,  denoted  by 
Cj,  j  =  1,2, ...,  J  h~^  =  .  Cj  contains  those  6(2j  +  l)^  — 12(2j  +  l)+8  =  24j^+2  grid 

points  Ynj  that  are  on  the  faces  of  the  cube  centered  at  ym  with  each  face  a  square  of  length 
2jh.  Then  we  have  jh  <  jy^^.  -  y^l  <  V3jh,ynj  G  Cj,  and  am,nj  <  {kjh)~°‘  = 

For  0  <  (5  <  1  arbitrary  close  to  0,  one  has  the  row  sum  estimate 

Eil  =  1  +  EU  <  1  +  E/=l(24j2  +  2)j-^-k-^-^ 


(70) 


(  1  +  fc3(l-5)-2a  <  p(l-5)-2a^  ^  <  3^  §«  <  1  -  5  <  1 


1  +  k  In  A:  <  1 , 


q;=§,  0<(5<1 


[  1  + <1,  a  >  §,  0  <  5  <  1. 

Also  this  is  true  for  any  point  ym  which  has  at  least  1/8  of  24j^  +  2  points  in  its  j-th 
cube  neighbor  Cj  for  j  =  1,  2, . . . ,  J  ~  k^~^ .  Hence  for  a  fixed  a  >  0  and  any  0  <  5  <  1 
arbitrary  close  to  0, 


(71) 

m=l  m=l  n=l 

Combining  (71)  with  (60),  we  have 

(l_g2)2pa^ 

(72) 


pi3{l-5)-a)^  a<| 


^3(1-5), 


a  > 


—  2 


Ni  > 

—k  ~ 


a  <  2 


Combine  (71)  with  (61),  we  have 


(73) 


(l-e2)2fc3(i-^),  a>i 


£-4A:2(3(i-5)-«),  a<| 


Nl  < 


a  >  I 


□ 


g-4^3(l-5)^ 

Replace  s5  by  5  in  each  of  the  above  cases,  we  complete  the  proof. 

Remark  3.1.  When  d  =  3,  if  either  dim(X)  =  dim(y)  =  3  or  X  and  Y  are  separated, 
there  are  0  <  c  <  c  <  oo  independent  of  k  such  that  c  <  ||G'(-,y)||2  <  c  uniformly 
in  y  G  Y.  For  the  same  set  of  Green’s  functions  G(x,  y^)  as  in  the  above  proof  and 
A  =<  G{-,ym),G{-,yn)  >,  we  have 


Nl 


cnl  <  E 


E  Am  >  (1  -  e^)  E 


m=l 


m=l 


m=l 


(74) 
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Hence  inequalities  (60)  is  replaced  by  the  following: 

<  m 

(75)  E  >  E  ^ 

m=l  m=l 


(1  -^cnl 

m 


[(1 


—  e 


)cn 


h]2 

kl 


m 


and  (61)  is  the  same.  The  estimate  of  —  X]m=i  Z]n=i  ®mn  is  amplified  by  (? 

from  the  previous  estimate  at  most.  So  the  same  results  are  true. 

When  d  =  2,  except  for  a  scaling  factor  k~'2  for  the  Green’s  function,  everything  else  is 
similar  to  the  ease  d  =  3.  By  the  definition  3.2  of  Nfj.,  it  is  independent  of  a  constant  scaling 
of  the  whole  set  of  vectors.  So  the  results  for  also  holds  for  A  =<  G{-,ym),G{-,yn)  >• 

Remark  3.2.  The  size  of  Y  can  be  scaled  for  the  result  of  Lemma  3.1.  If  a  >  ^,  trace 

estimates  (62),  (67),  and  (71)  gives  ^  (^^in  be  seen  easily  from  (60)  and 

(61)  that  ifY  is  scaled  to  oY ,  k  is  scaled  to  ak  in  those  estimates  (56),  (57)  for  Nl  and 
N respectively.  Otherwise,  the  only  factor  that  can  not  be  scaled  with  in  the  trace 
estimates  is  k~‘^°‘^  in  (62),  (66)  and  (70).  Hence,  in  addition  to  scale  k  to  ak,  there  is  an 
extra  factor  of  a?'^^  for  the  result  of  Lemma  3.1. 

Remark  3.3.  The  key  estimates  in  the  proof  of  Lemma  3.1  are  (60)  (or  (75) J,  (61),  and 
the  one  for  =  tr{A^A).  Although  the  estimate  of  tr{A^ A)  can  be  improved  by 

more  careful  estimates  of  each  row  sum  Yln=i  ^mn  taking  into  account  different  decorre¬ 
lation  rate  aceording  to  Theorem  2.1,  e.g.,  dividing  those  points  in  eaeh  j-th  square  (cube) 
neighbors  ofym  (see  Figure  6)  into  directional  cone  sections  according  to  whether  the  line 
eonnecting  y™  and  its  neighbors  intersecting  X  or  not,  which  gives  different  decorrelation 
rate  of  two  Green’s  functions  or  different  power  a  in  am,nj  ^  {kjh)~°‘  for  points  in  different 
cones,  as  long  as  there  is  a  solid  angle  cone  such  that  lines  connecting  those  neighbors  in 
that  section  and  ym  intersect  X ,  the  order  of  the  estimate  can  not  he  improved.  On  the 
other  hand,  whether  (60)  and  (61)  are  sharp  or  not  is  a  more  complicated  issue.  The  answer 
depends  on  the  variation  of  leading  eigenvalues  of  the  matrix  A  =<  G(-,  y^),  G(-,  yn)  >, 
which  depends  on  the  geometrie  setup  of  X  and  Y ,  and  the  choice  of  e. 

3.2.  Lower  Bound  and  Upper  Bound  Estimate  for  Approximate  Separability  of 
Green’s  Function.  Now  we  use  Lemma  3.1  to  prove  the  following  lower  bound  estimate 
for  approximate  separability  (14)  of  Green’s  function  for  Helmholtz  equation  (5)  in  the 
high  frequency  limit. 

Theorem  3.1.  Let  X,Y  be  two  compact  manifolds  embedded  in  R'^,d  =  2,3,  and  d  > 
dim(X)  >  dim(y)  =  s.  Assume  that  for  any  two  points  yi,y2  G  Y, 

I  <  G(-,yi),G'(-,y2)  >  I  <  ik\yi  -y2|)““ 


as  k\yi  -  y2|  oo 
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for  some  a  >  0.  If  there  are  fi{:x.)  £  L2{X),  gi{y)  £  L2{Y),  I  =  1,2, . . . ,  such  that 


(76) 


then 

(77) 


1=1 


<  e, 

L2{XxY) 


m> 


Cek‘^°',  a  <  |, 

,  a  >  §, 


for  any  0  <  5  <  1  and  arbitrary  close  to  0  as  k  ^  oo,  where  >  c(l  —  (Ce)^)^  for  some 
positive  constants  c  and  C  that  only  depend  on  X,  Y  and  n(x). 


Proof.  Without  loss  of  generality,  we  assume  y  is  a  line  of  unit  length,  a  unit  square  or 
a  unit  cube  for  s  =  1,2,3  respectively.  First,  put  down  a  uniform  coarse  grid  in  Y  with 
grid  size  h  =  k^~^,  0  <  (5  <  1  and  is  arbitrary  close  to  0.  The  grid  divides  Y  into  cells 
Ym,  m  =  1,2,...  Divide  each  coarse  cell  Ym  further  into  uniform  finer  cells 

of  size  Ji  <  k~^,  Ym,n,  n  =  1,  2, . . . ,  Njf  =  {h/hy.  See  Figure  7.  For  a  fixed  n,  Ym,n  is  in 
the  same  relative  location  in  each  coarse  cell  Ym-  The  center  and  volume  of  each  cell  17n,n 
is  ym,n  and  hf  respectively.  Define  Gh{x,y)  =  (^(x, y^.n), Vy  G  Ym,n  to  be  a  piecewise 
constant  function  in  y.  We  show  that  by  choosing  h<k~^  small  enough 


ih 


(78)  dy  |G(x,y)-Gfe(x,y)pdx  =  ^  ^  /  dy  |G(x, y)-G(x, ym,n)pdx  < 
JY  JX  —  JX 


If  X,  Y  are  disjoint,  assuming  n(x)  is  smooth  in  the  Helmholtz  equation  (5)  and  there  is  no 
caustics,  we  have  |VyG'(x, y)|  <  k  uniformly  in  x  G  X  and  y  £Y.  One  can  pick  h  <  k~^ 
small  enough  such  that 

|G'(x,y)  -  G(x,ym,„)|  <  e,  yx£X,y£Ym,n,  m  =  1,2, . . . ,  n]^  ,n  =  1,2, . . . ,  N^. 
and  get  (78). 


Figure  7.  Two  scale  decomposition  of  the  source  domain. 
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If  X  and  Y  are  not  disjoint  and  dim(X)  =  dini(y)  =  d  =  2  or  3,  G{-,y)  G  L2{X),\/y  G 
Y ,  we  show  that  there  is  still  h  <  k~^  small  enough  such  that 

(79)  [  |G'(x,y)  -  G'(x,ym,n)pdx  <  e^,  ^y  e  Ym,n, 

Jx 

which  implies  (78).  From  (12)  and  the  asymptotic  formulas  (7),  (8),  we  have 

(80)  [  |G(x,y)|2dx<r 

•iBr(y) 

with  the  constant  independent  of  k,  where  Br{y)  denotes  the  ball  with  radius  r  centered 
at  y.  Hence  there  are  balls  and  H.^(e)(ym,n)  with  radius  r(e)  ~  such  that 

(81)  f  |G'(x,y)pdx  <  e^,  [  |G(x, ym,n)|^dx  < 

dxnR,(,)(y)  dAnB,(,)(y,„,„) 

for  a  given  e,  Denote  X^  =  X  D  (H^(e)(y)  U  H.r(e)(ym,n))  and  Xf  =  X  \  X^.  For  y  G  Ym,n 
we  have 

fx  |G(x,y)  -  G{yL,ym,n)\^dx 

=  Jxc\G{xi,y)  -  G{x,ym,n)\‘^dx  +  J^JG{x,y)  -  G{x,y„i^n)\‘^dx  =  I  +  II 

Since  |VyG'(x,  y)|  <  max(/cr“^(e),  T“^(e)),  Vx  G  X^' ,  by  choosing  h  <  k~^  small  enough 
we  get  I  <  e^.  From  (81)  we  get  II  <  e^.  Hence  we  prove  (79). 

Let  the  linear  subspace  Sx  =  span{/p(x)}p^j^  C  L2{X),  then 

l|G(x,y)  -  Ps^G(x,y)||i^(^)dy  <  e^, 

where  Psx  is  the  projection  onto  Sx-  From  (78),  we  have 

J^UI  -  Psx)[G{^,y)  -  Gh{x,y)]\\l^^^-^dy  <  0{e^), 

where  I  is  the  identity  map,  from  (78)  we  have 

^  Iy  Wi^  -  Psx)Gh{^,y)]\\l^^x)dy 
~  Y^rn=l  Z]n=l/Y„  „  II  (7  ~ )^(^)  ym,n)  ||^2(X)'^y 
~  Sm=l  Yln=l  l|77(x,  ym,n)  —  PSxG{x.,  ym,n)||^2(X) 

=  {h/hyV  Y!^=1  YHU  l|G(x,ym,n)  -  PSxG[^Xm,n)\\l^(^x) 

~  Ylin=l  Ylni=l  ll^(^)  ym,n)  —  PSxGi^.,  ym,n)||^2(X) 


(82) 
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Assume 

(83) 


TV," 


\\G{x,y^^n)  -  PSxG{:^,ym,n)\\Yx)  =mm^  \\G{x,ym,n)  -  PSxGix,ym,n)\\Yx) 


m=l 

Then  there  is  a  constant  G  >  0  such  that 

Nf 


m=l 


(84) 


1 


h  \\G{^,ym,n)  -  PSxG{^,ym,n)\\Yx)  <  G' 

m=l 


2^2 


Since  ym,n  G  Y,m  =  1,2, .. .  forms  a  uniform  grid  with  grid  size  h  =  we 

can  apply  Lemma  3.1  to  get 


dim(S'x)  > 


c(l  -  (Ce)2)2A;2“  a  <  f 
c(l  -  {Ceffk^-^  a  >  f 


for  any  0  <  5  <  1  and  arbitrary  close  to  0  as  /c  — )•  oo,  where  C*  >  0,  c  >  0  are  some 
constants  that  depend  only  on  X  and  Y  and  n(x). 

□ 


Remark  3.4.  Theorem  3.1  presents  the  intrinsic  mathematieal  diffieulty  for  numerical 
computation  of  Helmholtz  equation  with  large  wave  number  k.  When  n  =  3,  the  same  lower 
bound  estimate  for  approximate  separability  holds  for  unnormalized  Green’s  function  too. 
Hence,  a  discretized  system,  such  as  the  discretized  kernel  in  boundary  integral  formulation 
or  inverse  of  the  matrix  corresponding  to  direct  discretization  of  the  PDE  and  its  off- 
diagonal  sub-matrices  corresponding  to  a  mesh  that  resolves  the  wavelength,  does  not  have 
low  rank  approximation  for  large  k  due  to  the  decorrelation  of  Green’s  funetions  caused  by 
fast  oscillations.  However,  when  n  =  2,  for  separated  X  and  Y,  |G(x,  y)|  <  /c“2  ^  V(x,  y)  G 
X  X  Y  which  approximate  separability  trivial  when  k  is  large  enough. 

Remark  3.5.  Although  it  seems  that  the  lower  bound  estimate  (77)  in  Theorem  3.1  shows 
a  weak  dependence  of  Nf  on  e  as  k  ^  oo,  the  choiee  of  e  can  affect  the  sharpness  of  the 
lower  bound  as  commented  in  Remark  3.3. 

Here  we  also  give  an  upper  bound  estimate  for  the  approximate  separability  of  Green’s 
function  in  the  high  frequency  limit.  The  intuition  is  that  the  Green’s  functions  with 
sources  located  within  a  wavelength  are  correlated.  So  use  the  linear  subspace  spanned  by 
Green’s  functions  sampled  uniformly  in  the  source  domain  with  a  grid  size  smaller  than 
wavelength  should  approximate  the  whole  family  of  Green’s  functions  good  enough. 

Theorem  3.2.  Let  X,Y  be  two  compact  manifolds  embedded  in  R'^,d  =  2,3  and  d  > 
dim(A)  >  dim(y)  =  s.  For  any  e  >  0  and  5  >  0,  there  are  fi(x)  G  L2{X),gi{y)  G 
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L2{Y),1  =  1,2, such  that 


(85) 


G(x,y)  -'^fi{x)gi{y) 


l=i 


<  e 

L2{XxY) 


as  k  ^  CO,  where  C  >  0  is  some  eonstant  that  depends  on  X,  Y  and  n(x). 


Proof.  Put  down  a  uniform  grid  in  Y,  which  is  assumed  to  be  a  line  of  unit  length,  a 
unit  square  or  a  unit  cube  for  s  =  1,2,3  respectively,  with  grid  size  h  =  k~^~^l‘^.  Denote 
the  grid  point  as  ym,in  =  1,2,  Denote  the  liner  subspace  Sx  = 

span{G(x,ym)},„ii  C  L2{X).  We  claim 

(86)  ||G'(x,y) -P5^G(x,y)||i2(A)  < 

for  k  large  enough,  where  |y|  denotes  the  volume  of  Y. 

If  X  and  Y  are  disjoint,  assuming  n(x)  is  smooth  in  the  Helmholtz  equation  (5)  and 
there  is  no  caustics,  we  have  |VyG(x,y)|  <  fc,  ||IlyG(x,  y)||  <  /c^,  where  the  bound  is 
uniform  in  X  and  Y.  Given  a  non-grid  point  y  £  Y,  G{x,  y)  can  be  approximated  by  a 
linear  interpolation,  which  is  convex  combination  of  Green’s  function  at  its  neighboring 
grid  points.  To  be  precise,  suppose  y  £  Y  lies  in  the  s  dimensional  simplex  with  vertices 
ymi ,  •  •  • ,  Yms+i  ■  Let  ry , . . . ,  ry"*"^  be  the  barycentric  coordinates  for  y,  i.e., 

S+1  5  +  1 

j=i  i=i 

Then  we  have  the  following  linear  interpolation 

s+l 

(87)  +(x,y)  -^YyG{x,ymj)\  <  \\DlG{x,y)\\h‘^  <  /c“^ 


and  hence  (86)  is  true  when  k  is  large  enough. 

If  X  and  Y  are  not  disjoint  and  dim(X)  =  dim(y)  =  d  =  2,3,  G{-,y)  £  L2{X).  For  a 
given  e,  there  are  balls  il.,-(^)(y)  and  H.,-(^)(ymj)  with  radius  0  <  r(e)  ~  centered  at  y 
and  ymj  ,j  =  1,  2, . . . ,  d  +  1  such  that 


(88) 


IB. 


■(o(y) 


|G'(x,y)|+x< 


8|y|(d  +  2)’ 


|G'(x,  y^Jpdx  < 


8|y|(d  +  2)' 


Alsowehave  |VyG'(x,y)|  <  max(/cr  ^{e),T  ^{e)),\\D^G(yL,y)\\  <  max{k^T  ^{e),kT  ^(e),r 

and  hence 

(89) 

5+1 

|G(x,y)  -  j;+G(x,y^J|  <  \\DlG{x,y)\\h^  <  max(A:2(i-^),  fci-27^-2(e),  fe-27^-3(e)). 
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for  any  x  G  X  \  -Br(e)(y)  U  (U^=i-BT(e)(ymj))-  By  decomposing  the  integration  in  (86) 
on  X  into  two  parts,  one  on  =  X  n  (-Br(e)(y)  U  (U^=i-B^(e)(ymj)))  and  the  other  on 
xf  =  X  \  (B^(,)(y)  U  (U^il5r(Q(ym,  ))),  we  have 

/x,  \G{^,y)  -  ES  ^y^(x,ym,)pdx 


Remark  3.6.  The  above  upper  bound  holds  for  unnormalized  Green’s  function  when  d  =  3. 

If  dim(X)  =  dim(y)  =  d  =  2,3,  since  the  Green’s  function  belongs  to  L2(X  x  Y), 
our  approximate  separability  estimates  in  L2  norm  is  valid  for  general  compact  domains 
X  and  Y,  disjoint  or  not.  However,  if  X  and  Y  are  disjoint,  which  is  the  case  for  most 
applications,  our  results  in  L2  norm  are  also  true  in  Loo(X  x  Y)  norm.  Since  Loo  norm  is 
stronger  than  L2  norm  in  a  compact  domain,  the  lower  bound  for  approximate  separability 
in  Theorem  3.1  immediately  extends  to  Loo  norm.  Also,  first  part  of  the  proof  in  Theorem 
3.2  directly  extends  to  L°°  norm.  We  summarize  these  two  results  below. 

Theorem  3.3.  Let  X,  Y  be  two  disjoint  compact  manifolds  embedded  in  R^,  d  =  2,3,  and 
d  >  dim(X)  >  dim(y)  =  s.  Assume  that  for  any  two  points  yi,y2  G  Y, 

I  <  G'(-,yi),G(-,y2)  >  |  <  (/c|yi  -y2|)"“  as  k\yi  -  y2|  ^  00 
for  some  a  >  0.  If  there  are  fi{:>L)  G  Loo(X),  gr;(y)  G  LaoiY),l  =  1,2,...,  Nf.  such  that 

y)-^fl  i^)9l  (y )  <  e,  Vx  G  X,  Vy  G  Y, 

1=1 


(93) 
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then 


(94) 


Cfk 


2a 


Cfk 


d—S 


a  <  f , 


a  >  §, 


for  any  0  <  5  <  1  and  arbitrary  close  to  0  as  k  ^  oo,  where  Ce  >  c(l  —  (Ce)^)^  for  some 
positive  constants  c  and  C  that  only  depend  on  X,  Y  and  n(x). 

Theorem  3.4.  Let  X,  Y  be  two  separated  compact  manifolds  embedded  in  d  =  2,3,  and 
d  >  dim(X)  >  dim(y)  =  s.  For  any  e  >  0  and  5  >  0,  there  are  fi{:>L)  G  L^{X),gi{y)  G 
LooiY),  I  =  1,2, ... ,  <  Ck'^~^^  such  that 


(95) 


as  k 


Nf 


G{^,y)  -  ^/i(x)5«(y) 


1=1 


<  e,  Vx  G  X,\fy  G  Y, 


oo,  where  C  >  0  is  some  constant  that  depends  on  X,  Y  and  n(x). 

Remark  3.7.  In  Theorem  3.2,  upper  bound  estimate  for  the  approximate  separability  of 
Green’s  function  for  the  Helmholtz  equation  in  the  high  frequency  limit  in  L2  norm  is 
derived  based  on  separable  approximation  using  linear  combination  (interpolation)  of  a  set 
of  Green’s  functions  with  sources  located  on  a  uniform  grid.  It  is  also  possible  to  obtain  an 
upper  bound  estimate  in  the  L2  norm  for  the  approximate  separability  of  Green’s  function 
for  the  Helmholtz  equation  in  the  high  frequency  limit  in  homogeneous  medium  for  an 
arbitrary  bounded  domain  based  on  separable  approximation  using  the  eigenfunctions  of 
Laplace  operator  and  Weyl’s  asymptotic  formula  [28]  for  the  eigenvalues. 

Suppose  G(x,  y)  is  the  Green’s  function  in  a  bounded  domain  Q  C  satisfying 

AxG(x,y)  +  A:2G(x,y)  =  (5(x  -  y),  x,yeQ, 

G(x,  y)  =  0,  X  G  dLl. 

Let  Um(x),  ||wm||L2(D)  =  1,  m  =  1,2, ...  be  the  normalized  eigenfunctions  for  the  Laplace 
operator 

(97)  Altm(x)  =  Ajtm(x),  X  G  O,  iim(x)  =  0,X  G 

with  eigenvalues  0  >  Ai  >  A2  >  ....  Hence  nm(x),m  =  1,2,...  are  also  the  normalized 
eigenfunctions  for  the  homogeneous  Helmholtz  operator  with  eigenvalues  Xm  +  k‘^,m  = 
1,2,....  Here  we  assume  the  domain  Ll  is  not  resonant,  i.e.,  Xm  +  kf  0,Vm.  Since 
Um{^)  forms  an  orthonormal  basis  for  L2{Ll)  and  G(x,y)  G  L2{Ll)  for  d  =  2,3,  one  has 
the  following  expansion 


(96) 


(98)  G{x,y)  ='^{Xm  +  k^)  ^Jtm(y)jtm(x). 

m=l 

The  Weyl’s  asymptotic  formula  gives 

(99)  I  A™  I 
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for  large  m,  where  |r2|  is  the  volume  ofQ..  Choose  a  large  enough  integer  M  >  ,  (5  >  0, 

then  \Xm  +  <  m~d  for  m  >  M .  For  any  e>  Q,  we  have 

„  M  oo 

(100)  /  /  |G(x,y)  -  ^  m~d<M^~d<e^ 

m=l  m=M+l 


when  k  is  large  enough  for  d  =  2,3.  Hence  one  can  use  eigenfunctions  of  Laplace  operators 
to  construct  a  separable  approximation  for  the  Green’s  function  of  Helmholtz  equation.  In 
particular  for  any  two  subdomains  X,Y  ofLl,  (100)  implies 


(101) 


G{-y^,y) +  ^«m(y)Hm(x) 

L2{XxY) 

<  G(x,y) -X;m=i(Am  +  A:2)“^Wm(y)Hm(x)  <e 


which  shows  that  <  k^~^^  for  any  5  >  0  in  the  high  frequency  limit. 


Remark  3.8.  It  can  be  seen  from  Theorem  3.1  and  3.2  as  well  as  Remark  3.7  that 
if  two  Green’s  functions  at  different  sources  decorrelate  fast,  <(?(•,  yi),  (?(•,  y2)  >  < 

(^|yi  “  y2|)~“  with  a  >  §,s  =  dim(y)  <  dim(X),  the  upper  and  lower  estimate  of  the 
approximate  separability  of  Green’s  function  for  Helmholtz  equation  is  sharp  in  the  high 
frequency  limit.  In  this  scenario,  a  set  of  Green’s  functions  with  sources  uniformly  dis¬ 
tributed  in  a  domain  or  a  set  of  leading  eigenfunctions  for  Laplace  operator  form  a  good 
basis  to  represent  an  arbitrary  Green’s  function  or  solution.  However,  the  representation 
is  not  much  compressible  and  hence  low  rank  approximation  in  numerical  computation  is 
not  feasible.  Moreover,  to  compute  a  set  of  densely  distributed  Green’s  functions  or  a  large 
number  of  eigenfunctions  of  Laplace  operator  of  order  at  least  0(/c®)  for  a  s-dimensional 
manifold  for  a  general  setup  in  practice  is  computationally  challenging.  In  certain  setups, 
explicitly  known  special  functions/basis  can  be  explored  to  develop  fast  algorithms  even  for 
dense  matrices  that  do  not  have  low  rank  approximation,  such  as  fast  Fourier  transform 
(FFT). 


4.  Examples 

In  this  section  we  apply  our  previous  general  results  to  a  few  setups  that  are  of  practical 
interests  in  3D.  Again  we  assume  that  X,  Y  are  two  compact  manifolds  embedded  in 
and  dim(A)  >  dim(y)  =  s.  First,  in  Section  4.1,  we  discuss  several  examples  for 
which  the  Green’s  function  for  Helmholtz  equation  is  not  highly  separable  in  the  high 
frequency  limit.  In  particular,  if  the  decorrelation  of  two  Green’s  functions  is  fast  enough, 
<  G(-,  yi),  G(-,  y2)  >  <  (^|yi  —  y2|)~“)  oi>  the  lower  bound  and  upper  bound  for  the 

approximate  separability  of  Green’s  function  for  Helmholtz  equations  in  the  high  frequency 
limit  are  sharp.  As  discussed  in  the  proof  of  Theorem  2.1,  Remark  2.3  and  Section  2.2,  if  a 
generic  ray  going  through  two  points,  yi,y2  G  Y ,  only  intersect  X  at  isolated  points,  then 
a  =  I  and  hence  both  lower  bound  and  upper  bound  are  sharp.  For  all  these  examples. 
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where  X  and  Y  are  given  and  fixed,  the  lower  bound  in  Theorem  3.1,  sharp  or  not,  implies 
that  there  is  no  low  rank  approximation  in  the  corresponding  discretized  system  due  to  the 
fast  decorrelation  of  Green’s  function  of  Helmholtz  equation  in  the  high  frequency  limit. 
However,  there  are  special  k  dependent  setups  for  X  and  V,  which  are  discussed  in  Section 
4.2,  where  the  Green’s  function  for  Helmholtz  equation  is  highly  separable  even  in  the  high 
frequency  limit.  High  separability  in  these  special  setups,  which  implies  availability  of  low 
rank  approximation  in  the  corresponding  discretized  system,  can  be  explored  to  develop 
fast  algorithms.  In  the  following  study,  all  constants  are  positive  and  only  depend  on  X,  Y 
and  n(x)  by  assuming  e  is  small  so  that  the  weak  dependence  of  those  estimates  in  Theorem 

3.1,  3.3  on  it  is  neglected. 

4.1.  Examples  of  not  highly  separable  Green’s  function. 

1)  X  and  Y  are  two  disjoint  compact  domains  in  dim(X)  =  dim(y)  =  s  =  3.  For  two 
points  yi,y2  G  Y,  we  can  only  claim  |  <  G(-,  yi),  G(-,  y2)  >x  |  <  (A:|yi  -  y2|)“^  in  general 
from  Theorem  2.1  since  for  any  point  y  GY  there  is  a  cone  with  y  as  the  vertex  and  a  solid 
angle.  A  segment  of  the  ray  (a  straight  line  in  homogeneous  medium)  connecting  y  and  a 
point  in  the  cone  stays  in  X  which  gives  a  ID  curve  of  stationary  phase.  Since  a  =  1  <  | 
our  lower  bound  and  upper  bound  estimates  are  not  sharp.  Theorem  3.1,  3.3  give  the  lower 
bound  while  Theorem  3.2,  3.4  give  the  upper  bound  <  k^~^^  for  any  h  >  0. 

2)  X  and  Y  are  two  disjoint  compact  surfaces  in  3D,  dim(A)  =  dim(y)  =  s  =  2.  This  is  a 
typical  scenario  for  boundary  integral  methods  when  X  and  Y  are  two  pieces  of  the  bound¬ 
ary  of  scatterers.  In  general,  the  ray  (straight  line  in  homogeneous  medium)  going  through 
two  points  yi,y2  G  Y  intersects  X  at  most  finite  number  of  times,  i.e.,  there  are  only 
isolated  stationary  points  for  the  oscillatory  surface  integral  of  <  (?(•,  yi),  (?(•,  y2)  >x  ,  so 

a  =  I  by  standard  stationary  phase  theory  (see  Remark  2.3).  From  Theorem  3.1-  3.4  we 
have  a  sharp  estimate 

Vh>0, 

as  k  ^  oo. 

Another  interesting  scenario  is  when  people  compute  the  direct  inverse  of  the  discretized 
linear  system  for  Helmholtz  equation  using  multi-frontal  method  in  which  the  full  linear 
system  is  reduced  to  smaller  but  dense  linear  systems  corresponding  to  unknowns  living 
on  planar  domains  such  as  those  depicted  in  Figure  8.  Low  rank  approximation  or  further 
skeletonization  for  these  smaller  but  dense  linear  systems  is  crucial  for  a  fast  numerical 
solver.  We  use  our  theory  to  give  the  approximate  separability  estimates  for  the  typical 
setups  in  the  high  frequency  limit.  Since  s  =  2,  the  upper  bound  is  <  k"^^^  for  any 
5  >  0.  Now  let’s  look  at  the  lower  bound.  We  first  look  at  three  typical  conhgurations 
in  homogeneous  medium  where  rays  are  (piecewise)  straight  lines:  (a)  if  X,  Y  are  two  dis¬ 
joint  coplanar  regions  as  shown  in  Figure  8(a),  the  least  rate  of  decorrelation  between  two 

Greens’s  function  is  a  =  he.,  <  G(-,  yi),  G(-,  y2)  >x  ^  (fc|yi  —  y2|)~2  by  Theorem  2.1 

when  a  segment  of  the  ray  going  through  yi,  y2  stay  in  A.  So  one  has  >  A:  as  A:  — )•  oo; 

(b)  if  X,  Y  are  two  disjoint  planar  regions  that  are  not  coplanar  nor  parallel  to  each  other. 
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e.g.  perpendicular  to  each  other  as  shown  in  Figure  8(b).  The  least  rate  of  decorrelation 
between  two  Greens’s  function  is  a  =  1  by  applying  the  standard  stationary  phase  result 
since  any  ray  going  through  two  points  yi,  y2  £  T  intersect  X  at  most  finite  times  (0  or  1). 

So  one  has  Nf,  >  >  0  as  /c  — 7-  oo;  (c)  if  X,  Y  are  two  planar  regions  in  parallel  as 

shown  in  Figure  8(c),  the  least  rate  of  decorrelation  between  two  Greens’s  function  is  a  =  | 

(in  general  since  any  ray  going  through  two  points  yi,y2  G  T  does  not  intersect  X)  or  1 
(if  part  of  the  boundary  dX  stays  on  a  level  set  of  the  phase  function  |x  —  yi|  —  |x  —  y2|), 
see  Remark  2.3.  So  one  has  V(5  >  0  as  /c  — )•  oo.  If  the  medium  is  inhomogeneous, 

a  ray  going  through  two  points  y i ,  y2  is  not  straight  line  and  intersect  a  planar  region  X 
only  a  finite  number  of  times.  So  a  >  1  and  we  have  the  sharp  low  bound  estimate 
0  in  general  for  the  above  three  configurations.  For  a  discretization  with 
a  fixed  ratio  of  grid  size  and  wavelength,  the  sharp  lower  bound  means  that  the  discretized 
matrix  and  its  sub  matrices  are  full  rank  modulo  a  constant. 


(a)  (b)  (c) 

Figure  8.  Two  planar  surfaces  in  3D. 


3)  X  G  is  a  compact  domain,  T  is  a  compact  smooth  curve  (dim(y)  =  s  =  1)  or  surface 
(dim(y)  =  s  =  2)  and  X,Y  are  disjoint.  Given  any  two  different  points  yi,y2  G  Y,  we 

have  <  G(-,  yi),  (?(•,  y2)  >x  <  (^|yi  “  y2|)~“  for  some  1  <  a  <  2  from  Theorem  2.1. 
Since  a  >  |,  from  Theorem  3.1-  3.4  one  has  the  following  sharp  estimates 

V<5>0, 

as  k  ^  oo. 

4)  X  is  a  2D  surface  or  ID  curve  and  T  is  a  ID  curve.  One  has  a  =  1  if  X  is  a  sur¬ 
face  or  a  =  2  if  -T  is  a  curve  using  standard  stationary  phase  theory  for  a  general  X  and 
two  points  yi,y2  G  Y.  So  for  both  scenarios  we  have  sharp  bounds 

V<5>0, 

as  k  ^  oo. 

4.2.  Highly  Separable  Cases.  Although  in  this  study  we  have  shown  that  the  lower 
bound  for  the  number  of  terms,  N^,  in  the  approximate  separability  (2)  of  the  Green’s 
function  G(x,  y)  for  Helmholtz  equation  in  the  high  frequency  limit  grows  with  certain 
power  of  k  (Theorem  3.1,  3.3)  in  general,  there  are  special  k  dependent  setups  of  X  and 
Y  where  G(x,  y)  is  highly  separable,  i.e.,  does  not  depend  on  k  and  depends  on  e 
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logarithmically.  These  special  setups  can  be  explored  to  develop  fast  algorithms  in  practice 
by  utilizing  low  rank  approximation. 

Let’s  assume  the  Green’s  function  is  of  or  can  be  approximated  by  the  form  G(x,y)  = 
y4(x,  where  yl(x,  y)  is  the  amplitude  function  and  (p{x,y)  is  the  phase  function. 

The  key  point  is  that  fast  decorrelation  or  almost  orthogothonality  between  Green’s  func¬ 
tions  with  different  source  locations  due  to  rapid  change  in  the  phase,  as  we  have  shown 
for  general  domains  X  and  Y,  are  not  present  in  those  special  setups.  Since  the  amplitude 
function  yl(x, y)  is  independent  of  k,  typically  one  needs  to  find  </>i(x)  and  02(y)  such 
that  k{(j){x,y)  —  4>i{x)  —  4>2{y))  is  uniformly  bounded  with  respect  to  x  £  X,y  G  Y  and 
k  and  has  a  highly  separable  approximation,  e.g.,  Taylor  expansion  in  x,y,  and  hence 
gifc(0(x,y)-(/>i(x)-</)2(y))  high  separable  approximation  due  to  the  fast  convergence  of 

Taylor  expansion  for  for  |2:|  <  C  <  oo.  Also  it  is  easy  to  see  that  the  phase  difference 
between  two  Green’s  function  at  different  sources  yi,y2  can  be  written  as 

yi)-(/>(x,  y2))  =  k{(l)2{yi)-4>2{y2))+k[{(l){x,  yi)-(/)i(x)-02(yi))-((/>(x,  y2)-</>i(x)-(/>2(y2))] 

which  is  a  constant  phase  shift  A:(i^2(yi)  —  4>2{y2))  plus  a  term  that  is  bounded  uniformly 
with  respect  to  x  G  A,  y  G  T  and  k  .  Hence  no  fast  oscillatory  integral  is  present  to 
decorrelate  two  Green’s  functions.  For  simplicity  3D  free  space  Green’s  function  (5)  are 
used  for  illustration  in  the  following  examples. 


1)  X,Y  are  two  line  segments  that  are  collinear  as  shown  in  Figure  9(a)  in  homogeneous 
medium 


(102)  <  G'o(-,yi),Go(-,y2)  >— 


||Go(-,yi)||2||Go(-,y2)||2 


Jk{y2-yi) 


lx  |x-yi||x-y2| 


rdx. 


There  is  no  fast  oscillatory  integral  to  decorrelate  two  Green’s  functions.  Denote  the  axis 
going  through  these  two  line  segments  as  r,  we  have 


(103) 


G'o(x,y) 


Q-ik(r^-ry) 

dvr  Ty  —  Tx 


In  this  trivial  case,  (j){x,y)  =  4>i{x)  +  4>2{y),  where  (/>i(x)  =  — rx,(/>2(y)  =  ry.  It  is  easy  to 
see  that  the  geometric  series  can  be  truncated  at  =  (log  log(47rpe)  to  get  an  a 

separable  approximation  for  any  e  >  0  independent  of  k,  where  lx  is  the  length  of  X  and 
p  is  the  distance  between  X  and  Y.  The  same  argument  can  be  applied  to  two  disjoint 
curve  segments  X  and  Y  that  lie  on  the  same  ray  in  heterogeneous  medium.  In  this  case 


(104)  <G(.,yi),G(.,y2)>= 


l|G(-,yi)||2||G(.,y2)|| 


.eiA:<^(yi,y2)  f  H(x,  yi)A(x,  y2)dx, 
!  Jx 


where  (/>(yi,y2)  is  the  travel  time  between  yi  and  y2  and  A(x,  yi),  A(x,  y2)  are  the  corre¬ 
sponding  amplitudes  in  geometric  optics  ansatz  as  discussed  in  Section  2.2. 

2)  X  and  Y  are  two  disjoint  thin  cylinders  around  a  line  as  shown  in  Figure  9(b).  This  3D 
setup  is  analogous  to  the  2D  setup  in  [21].  Numerical  test  will  be  presented  in  Section  5. 
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Figure  9.  Two  special  setups  of  X  and  Y  that  allow  highly  separable 
approximation  of  the  Green’s  function. 


Let  the  longitudinal  axis  be  r  and  the  other  two  orthogonal  coordinates  in  the  plane 
perpendicular  to  r  be  ^,77.  Let  x  =  {rx,Cx,'nx)  £  X  and  y  =  {ry,^y,riy)  £  Y.  Denote 
p  =  inf^ex,y&Yirx-ry)  and  r  =  supxgx,yeY  Assume  kr  <  =  ^  <  I-  Again 

in  this  case  =  —Vx,  </>2(y)  =  I’y  as  in  the  previous  case.  One  has 

/c|(/)(x,  y)  -  01  (x)  -  02(y)|  =  k{\x  -  y|  -  {ry  -  Vx))  <  2/cr  =  1 

Next,  we  give  an  explicit  separable  construction  using  asymptotic  expansions. 

(105) 

A:|x  -  y|  =  {rx  -  Vy)^  +  {^x  -  CyV  +  iVx  -  PyY 


=  Kry 


iyY  +  iVx 


—  -  (-l)-(2m)!  (fe 

y  n  —  2m)(m!)24™' 

m=l 


iy?  +  {rix  -  Pyf  T 
{ry  -  rxf^~^ 


Note  that  k^/ {^x  -  iyY  +  {Px  -  VyY  <  2/cr  <  1  and  So  second 

term  in  the  above  expression  can  be  bounded  by  a  geometric  series  S® 

(106) 


<  Go(-,yi),Go(-,y2)  >= 


1 


l|G'o(-,yi)l|2||G'o(-,y2)||2 


JKryi-ryY 


lx  |x-yi||x-y2| 


j-dx, 


where  0(x)  =  /c[(|x-yi|  -  {vy^  -  r^))  -  (|x-y2|  -  {vy^  -Ta;))]  and  |0(x)|  =  0(1),  Vx  e  X. 
Again  no  fast  oscillation  due  to  rapid  change  of  phase  is  present  in  the  integral.  Now  we 
show  explicitly  the  approximate  separation  based  on  the  expansion  (105). 

For  any  tolerance  e  >  0,  take  q  such  that  (^)^'^'''^(1  —  ^)~^  <  e-  So  only  the  hrst 
q  =  0(1  loge|)  terms  are  needed  in  the  summation  in  (105).  Denote  Q(x,y)  be  the  first  q 
term  expansion, 


g(x,  y)  =  k^(ix-iy?  +  {Vx-Vy?  E  (l-2i)(LT)2T 

m=l  '  /  \  / 


(-l)-(2m)!  ((C.  -  iy?  +  iVx  -  PyfT-^!'^ 


{ry  -  ra;)2™ 


-1 


Q{x,  y)  is  bounded  independent  of  k.  So 

g*fe|x-y|  _  giA:(ry-ra,)gi(Q(x,y)+e)  _  ^ik{ry-r^) ^iQ{yL,Y)  _|_ 
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Since  can  again  have  a.  p  =  0(|  loge|)  term  polynomial  expansion  for  any  tolerance  e  for 
a  bounded  x,  we  have 

1=0 


Let’s  look  at  each  term  in  the  summation.  Since  each  term  in  the  expansion  of  [(5(x,  y)]^ 
is  like 

■  / - 1^  +  {rix-riyf') 

k^JiCx  -  iy?  +  {Vx  -  vyY\  ^ ^ 

where  the  integer  m  ranges  from  I  to  {2q—l)l  and  m+l  is  even.  One  only  needs  to  keep  those 
m  <  2q  terms  in  the  expansion  because  =  0{e)  and  < 

1.  So  altogether  we  have  0{pq)  =  0(|  logep)  terms  of  the  form 

/-iQyN  ((^x  ~  ^y)'^  +  ivx  —  Vy)'^)^ 


where  0  <  /  <  m  <  |  loge|  are  integers.  Each  term 

((Cx  -  iy?  +  iVx  -  Vy?r  =  +  vl)  -  2Uy  -  2v.7Jy  +  ify  + 

can  be  expanded  into  =  0{\  logep)  separable  terms.  Further  more,  since  Vx  >  Vy  > 

0,  0(1  loge|)  leading  terms  are  needed  in  the  following  expansion  to  have  an  e  approxima¬ 
tion, 

i 


(108) 


{ry-r-r)  '  =  rj 


1-E 


oo  /  \  m 

r. 


m=l 


Hence  a  separable  e-approximation  of  requires  0(|  loge|®)  terms.  The  last  term  we 

need  to  make  a  separable  approximation  is  . 


|x-y| 


=  {fy  -  Tx) 


-1 


1  + 


(Cx  ^y)  {Vx  Vy) 
{tx  -  ryY 


2l 


^  ^  (-l)-(2(m  +  1))!  ((6  -  4)^  +  (r?x  -  r^yfT 

2(2m -b  l)((m -|- 1)!)24™'  {ry  —  rxY^ 

Both  summations  in  the  above  the  formula  can  be  truncated  at  0(|loge|)  terms  with  e 
error.  Each  term  in  the  summation  in  the  second  bracket  is  similar  to  the  term  in  (107) 
which  can  be  approximated  by  at  most  0(|  loge|^)  separable  terms.  So  can  also  be 

approximated  by  0(|loge|®)  separable  terms.  Combine  all  these  terms  together  we  have 
Nf.  <  0(1  loge|^^)  for  this  setup. 

3)  Here  we  include  two  setups  of  X  and  Y  that  has  have  been  proposed  and  used  for 
developing  fast  algorithms. 

(a)  In  [24,  6,  18]  fast  butterfly  algorithms  for  computing  highly  oscillatory  Fourier  integral 
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operators  and  boundary  integrals  for  Helmholtz  equation  were  developed.  The  key  idea 
is  a  dyadic  decomposition  of  two  domains  A,B  into  tree  structures  Ta,Tb,  from  root 
to  leaf,  and  a  recursive  paring  of  X  G  Ta  and  Y  ^  Tb  such  that  the  phase  function 
or  free  space  Green’s  function  restricted  on  X  x  Y  has  good  separability  property  or 
low  rank  approximation  in  discrete  setting.  The  condition  for  pairing  X  and  Y  is  that 
the  product  of  the  radius  of  X  and  Y  is  less  than  or  equal  to  For  this  setup,  the 
key  observation  is  that  one  can  construct  a  simple  separable  phase  function,  (/>(xo,y)  + 

(/>(x,  yo)  —  cf>(xo,yo)  that  can  approximate  the  original  phase  function  uniformly  well,  i.e., 
k\4’{x,y)  —  (j){xQ,y)  —  (/)(x,yo)  +  i^(xo,yo)|  is  uniformly  bounded  for  all  k,  x  £  X,y  G  Y , 
where  xo,yo  are  the  centers  of  X,Y  respectively.  Under  an  analytic  function  assumption 
for  </>(x,y),  k{(J){x,y)  -  (p{xo,y)  -  0(x,yo)  +  (/>(xo,yo))  can  be  approximated  by  a  Taylor 
expansion  with  0{  \  loge|)  terms  (see  the  proof  in  [6]). 

(b)  In  [8]  fast  directional  multilevel  algorithms  for  oscillatory  kernels  were  developed  based 
on  directional  low  rank  property  of  free  space  Green’s  function  of  Helmholtz  equation  on 
two  domains  Xr  and  that  satisfies  directional  parabolic  separation  condition:  is  a 

ball  of  radius  r  centered  at  a  point  c  and  X^  is  the  domain  containing  all  points  which 
are  at  a  distance  or  greater  from  c  and  belong  to  a  cone  centered  at  c  with  spanning 
angle  Here  r  can  be  thought  in  unit  of  wavelength  A  =  ^.  For  this  setup,  the  phase 
function  |x  — y|  =  ^|Ax  — Ay|  can  be  uniformly  approximated  by  x  •  (x  —  y)  =  |x|— x-y  = 

|x|  —  (x  —  ^)-y  —  Gy  which  can  be  further  approximated  by  |x|  —  1  •  y  uniformly,  where 
X  =  and  I  is  the  direction  the  cone  centered  at.  In  another  word,  let  (p{x,y)  =  |x  —  y| 

and  <?i>i(x)  =  |x|,(/)2(y)  =  —I  -  y,  then  |(/)(x,y)  —  (/>i(x)  —  4>2{y)\  is  uniformly  bounded  and 
has  an  0(log  |e|)  term  separable  approximation  with  any  e  >  0  error. 

4)  Here  we  present  a  scaling  argument  to  show  a  high  separable  asymptotic  regime 
for  two  general  disjoint  domains,  X  and  Y,  as  k  ^  oo.  Then  we  use  this  asymptotic 
argument  to  show  that  the  condition  for  the  setup  of  butterfly  algorithm  is  tight.  Let 
r{X),r{Y)  denote  the  diameters  of  X^Y  respectively.  Denote  the  distance  between  X 
and  Y  to  be  dist{X,Y).  Without  loss  of  generality,  let  r{Y)  <  r{X)  and  r{X)  =  0(1). 
Actually  the  size  of  X  is  not  restricted  since  it  can  always  be  scaled  to  0(1)  by  scaling 
X  to  for  the  Helmholtz  equation  (5).  Assume  r{Y)  <C  dist{X,Y).  From  the  scaling 
argument  in  Remark  2.2,  the  rate  of  decorrelation  of  two  Green’s  functions  is  rescaled 

to  <  G(-,yi),G(-,y2)  >  <  ,Vyi,y2  G  Y  for  some  a  >  0.  So  if  r(Y)  < 

^ dist{x^ ^  there  is  no  fast  decorrelation  of  two  Green’s  function.  On  the  other  hand,  if 

r(y)  >  \/^  >  0,  one  can  put  a  grid  with  a  grid  size  a  little  larger  than 

in  Y  so  that  they  become  more  and  more  decorrelated  as  A:  — )■  oo.  At  the  same  time,  the 
number  of  grid  points  in  Y  becomes  larger  and  larger  as  /c  — )■  oo.  Use  a  similar  argument 
as  in  Lemma  3.1  one  can  show  a  lower  bound  estimate  for  approximate  separability  which 
grows  with  some  power  of  k.  The  above  discussion  can  be  further  scaled  with  dist{X,  Y). 

For  example,  if  dist{X,Y)  =  0{k),  the  size  of  Y  can  be  0(1). 
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In  a  typical  setup  for  butterfly  algorithm  [24,  6,  18],  there  are  two  domains  A,  B  whose 
sizes  are  of  0(1)  and  dist{A,B)  =  0(1).  For  example,  A,B  may  be  disjoint  boundaries  of 
scatterers  in  boundary  integral  methods.  Dyadic  decomposition  of  A,  B  gives  tree  struc¬ 
tures  Ta,Tb  with  L  =  0{logk)  levels,  where  the  roots  are  A,B  and  the  leaf  nodes  at 
level  L  are  of  size  0{k~^).  Then  the  interaction  between  A  and  B  through  a  highly  os¬ 
cillatory  kernel  function,  e.g.,  the  Green’s  function  for  Helmholtz  equation  or  a  Fourier 
integral  operator,  is  computed  based  on  a  recursive  pairing  of  nodes  X,  A  D  X  £  Ta, 
and  Y ,  B  D  Y  £  Tb,  such  that  the  level  of  X,  1{X),  and  the  level  of  Y,  1{Y),  satisfy 
1{X)  + 1  (Y)  =  L.  The  key  observation  is  that  the  interaction  between  X  and  Y  has  a  low 
rank  approximation.  It  was  shown  in  [6]  that  the  low  rank  approximation  is  guaranteed  if 
r{X)r{Y)  =  0{k~^)  which  is  implied  by  condition  1{X)  +  1{Y)  =  L  (plus  some  analyticity 
condition  on  the  phase  function).  It  can  be  easily  seen  that  high  separability  for  the  setup 
for  butterfly  algorithms  falls  into  the  asymptotic  regime  discussed  above.  The  requirement 
of  analyticity  of  the  phase  function  is  equivalent  to  requiring  dist{X,  Y)  >  c  >  0  for  some 
fixed  c  and  the  condition  r{X)r{Y)  =  0{k~^)  implies  the  smaller  domain  Y  (or  X)  satisfies 
r(y)(or  r{X))<  0{k~2).  In  particular,  this  condition  is  barely  satisfied  when  r(X)  r{Y) 
or  1{X)  =  1{Y).  If  r{A)r{B)  =  0{k^)  for  any  5  >  0  or  dist{X,Y)  — ^  0  as  A:  — >  oo,  then  the 
condition  is  violated  and  low  rank  approximation  is  not  valid  as  discussed  above. 

4.3.  Approximate  separability  with  boundary  condition.  So  far  we  have  shown 
approximate  separability  estimates  for  Green’s  function  of  Helmholtz  equation  in  high 
frequency  limit  either  in  the  whole  space.  Here  we  present  an  example  when  boundary 
and  reflection  are  present.  First  we  study  approximate  separability  of  Green’s  function 
for  Helmholtz  equation  in  half  space  with  homogeneous  Dirichlet  boundary  condition  as 
illustrated  in  Figure  10(a).  The  Green’s  function  for  the  upper  space,  denoted  by  Gi(x,  y), 
can  be  explicitly  constructed  from  free  space  Green’s  function,  denoted  by  G'o(x,  y), 

Gi(x,y)  =  G'o(x,y)  -  Go(x,y), 

where  y  is  the  mirror  image  of  y  with  respect  to  the  boundary.  Decorrelation  of  Gi  can 
be  deduced  from  that  of  Gq.  Given  two  disjoint  compact  domains  X  and  Y  in  the  upper 
half  space  and  two  points  yi,y2  G  Y,  if  the  line  lyl  connecting  these  two  points  intersects 
with  X,  or  none  of  the  the  lines  iZ? ,  iZZ ,  intersects  with  X,  we  have 

7  ^  i  7  Yl  7  Yi  7 

(109) 

l<  G'i(-,yi),G(-,y2)  >| 

=|<Go  (x,y  i)  ,Go  (x,y2)>-<Go  (x,y  i)  ,Go  (x,y2)>-<Go  (x,y  j)  ,Go  {x,y2)>+<Go  (x,y  j)  ,Go  (x,y2)>l 
<  l<  G'o(-,yi),Go(-,y2)  >| 

since  |yi  -y2|  <  min(|yi  -y2l>  iFi  “  y2Myi  - y2l)  and  dist{yj,X)  <  dist{yj,X),j  =  1,2. 
Otherwise,  decorrelation  of  Gi(x,  yi)  and  Gi(x,  y2)  may  be  slower  than  that  of  Go(x,  yi) 
and  Go(x,  y2).  In  general,  for  two  disjoint  compact  domains  X  and  Y,  there  is  a  cone 
with  solid  angle  at  each  point  y  £  Y  with  y  as  the  vertex  such  that  a  line  connecting  y 
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and  any  other  point  in  the  cone  will  intersect  X.  Correlations  between  Green’s  functions 
of  these  points  are  the  leading  terms  in  the  estimate  of  tr{A^ A)  as  discussed  in  Remark 
3.3.  Hence  the  relation  (109)  implies  that  the  proof  and  result  in  Lemma  3.1  and  the 
lower  bound  estimates  in  Theorem  3.1  and  Theorem  3.3  for  X  and  Y  without  boundary 
hold  for  this  case  too.  It  can  also  be  easily  seen  that  separable  approximation  of  Green’s 
function  Gi(x, y)  for  x  £  X,y  £  Y  can  be  obtained  by  combining  separable  approximation 
of  Cro(x, y)  for  x  £  X,y  £  Y  and  that  of  G'o(x, y)  for  x  £  X,y  £  Y ,  where  Y  is  the  mirror 
image  of  Y  with  respect  to  the  boundary.  So  the  upper  bound  estimates  in  Theorem  3.2 
and  Theorem  3.4  hold  too  here.  Although  the  asymptotic  formulas  of  the  lower  bound  and 
upper  bound  estimate  for  approximate  separability  of  Green’s  function  in  this  case  is  similar 
to  those  in  free  space,  the  number  of  terms,  N^,  needed  for  a  separable  approximation  for 
the  Green’s  function  with  a  given  error  tolerance  e  >  0  defined  in  (2)  will  increase  due 
to  reflection  at  the  boundary  in  general.  It  can  be  argued  by  the  following.  The  Green’s 
function  Gi(x,y)  viewed  as  a  family  function  in  X  parametrized  by  y  G  T  is  composed 
of  a  family  of  free  space  Green’s  function  Go(x,y)  with  y  in  a  larger  domain  Y  UY.  Or 
geometrically,  instead  of  there  is  one  ray  connecting  any  two  points  x  £  X  and  y  G  T  in 
free  space,  now  there  are  two  rays  with  the  presence  of  a  reflection  boundary.  The  phase 
function  becomes  more  complicated  and  needs  more  terms  in  a  separable  approximation. 
As  the  distance  from  the  boundary  to  X  and  Y  becomes  larger  and  larger,  the  effect  from 
the  boundary  becomes  less  and  less.  In  general,  one  can  expect  that  presence  of  reflection 
and  scattering  can  cause  the  increase  of  N^,  the  number  of  terms  needed  for  a  separable 
approximation  for  the  Green’s  function  defined  in  (2). 


Figure  10.  Half  space  setups  with  a  reflection  boundary. 

Here,  we  give  a  more  explicit  study  of  the  boundary  effect  on  a  particular  setup  as 
depicted  in  Figure  10(b),  where  X  and  Y  are  two  disjoint  collinear  line  segments  parallel  to 
the  boundary.  Without  the  reflection  boundary,  Go(x,  y),  x  G  A,  y  G  F  is  highly  separable 
which  is  discussed  as  the  first  special  case  in  Section  4.2.  However,  if  the  distance  from  the 
boundary  to  X,  Y  is  small  compared  to  the  wavelength  and  the  separation  distance  between 
X  and  Y,  the  special  setup  of  two  thin  cylinders  in  Section  4.2  or  study  in  [21,  9]  shows 
that  Gi(x,  y),x  G  X,y  £  Y  is  still  highly  separable.  This  is  the  key  observation  for  the 
low  rank  approximation  used  in  the  sweeping  preconditioner  for  Helmholtz  equation  in  2D 
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in  [9].  However,  this  highly  separable  property  does  not  hold  if  the  boundary  is  not  close 
enough.  Also  high  separability  does  not  hold  in  an  analogous  setup  in  3D.  As  discussed 
in  case  2)  in  Section  4.1,  even  for  two  general  disjoint  compact  co-planar  2D  regions  X,Y 
without  the  boundary  the  number  of  terms  needed  for  a  separable  approximation  is  at  least 
0{k)  as  /c  — )■  oo. 

We  can  extend  the  above  study  further  to  Helmholtz  equation  in  a  waveguide,  see  Figure 
11.  The  Green’s  function  in  the  waveguide  with  a  source  at  y  is  an  infinite  sum  of  homoge¬ 
nous  Green’s  functions  with  sources  that  are  mirror  images  of  each  other  with  respect  the 
top  and  bottom  boundary,  i.e., 

CXD  OO 

(110)  G2(x,  y)  =  Go(x,  y)  +  ^  (-l)”^Go(x,  y+)  ^  (-l)”^Go(x,  y“ ), 

m=l  m=l 

where  yf  is  the  mirror  image  of  y  with  respect  to  the  top  boundary,  y^,  m  =  2, 3, ...  is  the 
mirror  image  of  y^-i  with  respect  to  the  top  boundary.  Similarly,  yj“  is  the  mirror  image 
of  y  with  respect  to  the  bottom  boundary  and  y“  ,  m  =  2,  3, ...  is  the  mirror  image  of  y^-i 
with  respect  to  the  bottom  boundary.  Due  to  the  two  reflection  boundaries,  the  Green’s 
function  G2(x,y),y  G  T  in  the  waveguide  involves  free  space  Green’s  function  Go(x,  y) 
with  y  belonging  an  infinite  union  of  Y  and  its  images  and  their  images  with  respect  to 
the  top  and  bottom  boundary  respectively.  In  another  word,  there  is  one  direct  ray  in  free 
space  and  inhnite  many  broken  rays  due  to  the  waveguide  setup  connecting  any  two  points 
X  G  A  and  y  G  T.  Of  course  the  decay  factor  will  play  a  role  here  too.  For  the  setup 
of  two  collinear  line  segments  X  and  Y  shown  in  Figure  11(b),  highly  separability  does  not 
hold  even  if  both  boundaries  are  very  close. 

reflection  boundary 

reflection  boundary 

X  Y 


y 

reflection  boundary  reflection  boundary 

(a)  (b) 

Figure  11.  Helmholtz  equation  in  a  waveguide. 


5.  Numerical  tests 

Here  we  show  a  few  numerical  tests  related  to  our  analysis  in  previous  sections.  In  all 
our  numerical  tests  free  space  Green’s  functions  are  used.  Our  computational  grid  size  h 
resolves  the  wavelength  A  =  2n/k,  h  =  A/15  in  2D  and  h  =  A/13  in  3D. 
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Example  1  We  test  the  decorrelation  of  two  Green’s  functions,  <  Gi(-,  yi),  G'i(-,  y2)  >x 
where  X  is  a  compact  domain.  In  this  test,  we  show  results  for  k  starting  from  50  with  an 
increment  of  5. 

Figure  12  shows  the  decorrelation  of  two  Green’s  functions  in  2D.  Here  the  domain  X 
is  a  disc  centered  at  (0,0)  with  radius  0.4.  In  Figure  12(a)  the  two  points  yi,y2  he  on 
x-axis.  So  the  line  through  yi,  y2  intersects  X  and  hence  there  is  a  line  of  stationary  points 
for  <  Gi(-,  yi),  Gi(-,  y2)  >x-  In  Figure  12(b),  the  line  through  yi,y2  is  parallel  to  y-axis 
which  does  not  intersect  X.  Hence  there  is  no  stationary  point  for  <  G'i(-,  yi),  G'i(-,  y2)  >. 
As  shown  in  the  proof  of  Theorem  2.1,  two  Green’s  functions  in  this  case  decorrelate  much 
faster  than  the  afore  mentioned  case. 

Also  two  Green’s  functions  decorrelate  faster  when  |yi  —  y2|  becomes  larger  in  general 
as  shown  in  Theorem  2.1,.  It  can  be  seen  that  as  yi,y2  are  further  away  from  A,  two 
Green’s  functions  become  more  correlated  due  to  scaling  of  the  gradient  and  the  Hessian 
of  the  phase  function  in  term  of  the  distance  from  yi,  y2  to  X  as  explained  in  Remark  2.2. 


Figure  12.  Decorrelation  of  two  Green’s  functions  in  2D. 


Figure  13  shows  a  test  with  homogeneous  Dirichlet  boundary  condition.  The  boundary 
is  an  infinite  line  that  is  parallel  to  x-axis  located  at  y  =  —d.  In  Figure  13(a),  X  is  again 
a  disc  as  before.  One  can  see  that  the  boundary  condition  does  affect  the  correlation 
between  two  Green’s  functions.  However,  the  asymptotic  behavior  is  similar  to  the  one 
without  boundary  condition  as  A:  — )•  oo.  In  Figure  13(b),  A  is  a  line  segment  on  x-axis 
between  [—0.5, 0.5],  which  is  co-linear  with  points  yi,y2-  With  no  boundary  conditions, 
the  two  Green’s  functions  are  just  different  by  a  constant  phase  and  hence  fully  correlated. 
When  the  boundary  is  present  and  is  not  too  close  or  too  far  away  from  the  two  points 
and  A,  we  do  see  more  decorrelations  due  to  the  boundary  boundary.  When  the  boundary 
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is  either  very  close  or  the  boundary  is  very  far,  the  two  Green’s  functions  become  highly 
correlated  as  explained  in  Section  4.3. 


x9 
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BC  atd 
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d=5,  "x" 
d=10, 


'  •  •  .  9000000065^  , 


O  000^0 
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®  d=0.5,  o 

d=5,  "x" 
d=10, 


150 


200 


(a) 


(b) 


Figure  13.  Decorrelation  of  two  Green’s  functions  in  2D  with  a  boundary. 

Figure  14  shows  corresponding  tests  in  3D.  The  behavior  in  3D  is  similar  to  those  in 
2D.  The  domain  X  is  a  ball  centered  at  the  origin  with  radius  0.4.  Figure  14(a)  shows  a 
test  where  the  two  points  yi,y2  lie  on  x-axis,  hence  the  line  through  yi,y2  intersects  X 
and  two  Green’s  functions  decorrelate  relatively  slow.  Figure  14(b)  shows  the  effect  of  a 
Dirichlet  boundary  condition. 


Figure  14.  Decorrelation  of  two  Green’s  functions  in  3D. 

Example  2  Here  we  present  singular  value  decomposition  (SVD)  pattern  for  matrix 
G(xj,yj),  where  Xj,yj  are  uniformly  distributed  points  in  X  and  Y  respectively  with  a 
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grid  size  that  resolves  the  wavelength.  3D  free  space  Green’s  function  are  used  in  the  test. 
In  these  tests,  we  show  results  for  k  ranging  from  10  to  150  with  an  increment  of  5.  In  the 
following  figures,  solid  lines  show  the  number  of  leading  singular  values  that  is  larger  than 
e  vs.  the  wave  number  k  while  dotted  lines  shows  ^  defined  by  Dehnition  3.2.  In  all  tests 
we  use  three  threshold  values  e  =  10“^,  10“^,  10“®. 

Figure  15(a)  is  the  result  for  X  =  {{x,y,z)\x  G  [— 0.5,0.5],y  =  0,2;  =  0}  and  Y  = 
{(x,  y,  z)\x  =  =  0.2, 2;  G  [—0.5,  0.5]}.  Figure  15(b)  is  the  result  for  X,  a  sphere  centered 

at  origin  with  radius  0.5,  and  Y  =  {{x,y,z)\x  G  [0.6, 1.6], y  =  0.6,2;  =  0}.  Linear  growth 
as  predicted  by  our  analysis  is  clearly  seen  for  both  cases. 


(a) 


(b) 


Figure  15.  Leading  singular  values  vs. 


wave  number  for  dim(y)  =  1. 


Figure  16  gives  an  example  of  two  square  domains  X,  Y  of  length  0.4  for  each  side  with 
different  relative  position  corresponding  to  the  three  setups  as  discussed  in  Section  4.1  case 
2)  and  demonstrated  in  Figure  8.  In  Figure  16(a),  X,Y  are  side  by  side  and  coplanar 
with  minimum  distance  0.1  between  them.  In  Figure  16(b),  X,Y  are  side  by  side  and 
orthogonal  to  each  other  with  minimum  distance  0.1  between  them.  In  Figure  16(a),  X,Y 
are  parallel  with  distance  0.3  between  them.  As  analyzed  in  Section  4.1,  case  (a)  has 
the  slowest  rate  of  decorrelation  between  two  Green’s  function  and  hence  also  the  slowest 
growth  of  the  number  of  leading  singular  values  and  while  case  (c)  has  the  fastest  rate 
of  decorrelation  between  two  Green’s  function  and  hence  also  the  fastest  growth  of  the 
number  of  leading  singular  values  and  among  the  three  cases  as  A:  — >•  00.  At  least  linear 
growth  is  observed  in  case  (a).  In  both  case  (b)  and  (c)  a  quadratic  growth,  as  predicted 
by  sharp  lower  bound  and  upper  bound  estimates  and  analyzed  in  Section  4.1  case  2,  is 
observed. 

Figure  17(a)  shows  an  example  of  two  spheres  of  radius  0.2  with  a  separation  distance 
of  0.2  between  the  two  spheres.  At  the  maximum  wave  number  k  =  150,  20,000  points 
are  randomly  distributed  on  the  surface  of  each  sphere  with  the  number  of  points  being 
proportional  to  wave  number  for  wave  numbers  in  between.  Again  one  sees  that  the  number 
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(a) 


(b) 


(c) 


Figure  16.  Leading  singular  values  vs.  wave  number  for  two  squares  in  3D. 


of  leading  singular  values  grows  quadratically  as  analyzed  in  Section  4.1  case  2).  Figure 
17(b)  shows  an  example  of  two  thin  cylinders  as  illustrated  in  Figure  9(b).  The  radius 
of  each  cylinder  is  half  wavelength  and  the  length  of  each  cylinder  is  0.4.  The  separation 
distance  between  the  two  cylinders  is  one  wavelength.  The  SVD  pattern  agrees  with  our 
high  separability  result  in  Section  4.2  nicely. 

6.  Conclusion 

In  this  work,  approximate  separability  of  Green’s  functions  of  Helmholtz  equations  in 
the  high  frequency  limit,  which  has  direct  implication  for  low  rank  approximation  for 
the  corresponding  discretized  linear  system,  is  studied  in  details.  By  characterizing  the 
decorrelation  rate  of  two  Green’s  function  due  to  fast  oscillations  in  various  situations  and 
showing  a  tight  dimension  estimate  for  the  approximation  of  a  set  of  almost  orthogonal 
vectors,  we  prove  an  explicit  and  sharp  asymptotic  formula  for  the  lower  bound  for  the 
number  of  terms  needed  for  a  separable  approximation  of  Green’s  function  in  the  high 
frequency  limit.  It  gives  a  rigorous  mathematical  argument  for  the  complexity  and  difficulty 
for  solving  high  frequency  Helmholtz  equation  numerically.  Application  to  setups  that 
are  commonly  used  in  practice  is  presented.  Numerical  tests  show  full  agreement  with 
the  analysis.  Development  of  efficient  numerical  algorithms  based  on  this  study  will  be 
investigated  in  the  future. 
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(a)  (b) 

Figure  17.  SVD  pattern  for  two  spheres  and  two  thin  cylinders. 
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